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ABSTRACT 


This  research  is  directed  at  investigation  of  regional  wave  excitation  and  propagation  in 
complicated  crustal  waveguide  using  numerical  simulations.  The  report  covers  three  major 
approaches:  a  finite-difference  plus  slowness  analysis  method  for  investigating  near-source  Lg- 
wave  excitation  and  partitioning;  the  half-space  screen  method  for  simulating  high  frequency  Lg- 
wave  propagation;  and  a  boundary  element  method  for  calculating  Lg-wave  propagation  in  a 
crustal  waveguide  with  surface  topography. 

(1)  A  finite-difference  modeling  plus  slowness  analysis  method  is  developed  to  investigate  the 
near  source  Lg-wave  excitation  and  energy  partitioning.  The  method  allows  a  very  fine  structure 
to  be  used  for  simulating  the  near-source  processes.  The  slowness  analysis  is  used  for  tracking 
where  energy  will  be  partitioned  into  the  long-range  propagation  regime.  The  high  efficiency  of 
the  method  allows  the  examination  of  many  source-model  combinations.  The  near  source  P-to-S 
conversion  in  the  presence  of  small-scale  random  heterogeneities  is  tested  as  a  mechanism  for 
generating  the  near  source  S-wave  energy  from  an  explosive  source.  The  numerical  results  reveal 
that  the  depth  of  the  source  and  the  depth  of  the  scattering  region  have  strong  effects  on  the  P-to- 
S  conversion  and  the  partitioning  of  energy  into  trapped  and  leaking  signal.  The  excitation 
spectrum  for  the  trapped  energy  from  the  S*-wave  is  also  investigated.  The  modeling  shows  that 
the  excitation  is  generally  favorable  for  lower  frequencies  and  shallow  source  depth.  We 
simulate  the  amplitude  ratio  between  the  P-  and  trapped  S-waves  for  a  relatively  broad  frequency 
range  using  velocity  models  with  and  without  lateral  velocity  variations.  For  a  typical  explosion 
source  and  a  shallow  earthquake  source  the  results  show  a  general  behavior  similar  to  that  of  the 
observations. 

(2)  The  half-space  screen  method  has  been  further  improved  to  model  Lg-wave  propagation  in 
the  complex  crustal  waveguides,  including  irregular  topography,  large-scale  geological  structures 
and  small-scale  random  heterogeneities.  The  screen  method  has  been  extended  to  the  P-SV  case. 
Both  P-  and  S-wave  quality  factors  are  incorporated  into  P-SV  and  SH-wave  propagators  to 
handle  the  effects  of  intrinsic  attenuation  on  regional  Lg.  The  resulting  propagator  provides  a 
useful  tool  for  investigating  the  relationship  between  Lg  wave  attenuation  and  crustal  waveguide 
properties.  Different  attenuation  mechanisms  including  intrinsic  attenuation,  leakage  due  to 
rough  topography  and  Moho  discontinuity,  and  scattering  attenuation  due  to  small-scale  random 
heterogeneities  are  tested  and  investigated  separately.  Then,  their  combined  effects  on  Lg-wave 
propagation  in  near-realistic  crustal  waveguides  are  estimated  and  the  corresponding  apparent  Lg 
quality  factors  are  calculated.  The  synthetic  apparent  Lg  Q  and  its  frequency  dependence  fall  into 
the  range  of  actual  observations. 

(3)  An  indirect  boundary  element  method  is  developed  to  calculate  the  two-dimensional  P-SV 
elastic  response  for  crustal  waveguide  model.  The  method  allows  us  to  accurately  handle  the 
wave  propagation  in  models  with  irregular  topography.  The  validity  and  accuracy  of  the  method 
are  checked  by  comparing  the  boundary  element  calculations  with  the  results  from  other 
numerical  methods.  Numerical  simulations  with  this  scheme  show  that  rough  topography  can 
scatter  the  P  and  Rayleigh  wave  and  attenuate  the  energy  propagation  in  the  waveguide.  In  order 
to  simulate  long-range  high-frequency  propagation  of  regional  waves  and  to  avoid  calculating 
huge  matrixes,  a  connection  technique  is  proposed  and  tested.  Using  the  technique,  a  long 
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waveguide  can  be  divided  into  relatively  short  sections,  and  the  boundary  element  method  can  be 
used  section  by  section  to  calculate  the  effect  of  rough  topography  on  wave  propagation  at 
extended  regional  distances. 
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1.  RESEARCH  OBJECTIVES 


With  the  current  emphasis  on  global  monitoring  for  low-yield  nuclear  tests,  regional  seismic 
phases  such  as  Lg  and  Lg  coda  have  become  very  important  for  magnitude  and  yield  estimation 
of  underground  nuclear  tests,  (e  g.,  Nuttli,  1986;  Xie,  et  al.,  1996;  Patton,  2001).  In  addition, 
various  P/S-type  amplitude  ratios  for  high  frequency  regional  phases  (e  g.,  Pn/Sn,  Pn/Lg,  Pg/Lg, 
Pg/Sn)  have  become  important  for  event  discrimination  (e  g.,  Taylor  et  al.,  1989;  Kim  et  al., 
1993,  1997;  Walter  et  al.,  1995;  Fisk  et  al.,  1996;  Taylor,  1996;  Taylor  and  Hartse,  1997;  Hartse 
et  al.,  1997;  Fan  and  Lay,  1998a-c;  Xie  and  Patton,  1999).  The  applications  of  regional  phases 
for  yield  estimation  and  event  discrimination  are  largely  based  on  empirical  approaches,  and 
while  very  promising  in  many  cases,  there  is  great  need  for  quantification  of  how  and  why  they 
work,  and  equally  importantly,  how  and  why  they  may  fail  under  certain  circumstances.  These 
questions  are  closely  related  to  the  Lg  wave  excitation  and  propagation.  Due  to  the  complex 
processes  involved,  it  is  difficult  to  separate  the  contribution  of  individual  mechnisms  from  the 
observed  data  empirically.  Numerical  modeling  approaches  are  thus  of  great  importance  for 
investigating  the  excitation  and  propagation  of  regional  phases. 

A  major  goal  of  the  UCSC  Regional  Wave  Synthetic  Seismogram  research  program  is  to  develop 
computationally  viable  techniques  for  investigating  the  excitation  and  path  effects  for  high- 
frequency  regional  waves.  The  techniques  that  have  been  developed  include:  ( 1 )  The  finite- 
difference  simulation  combined  with  slowness  analysis  method  (Xie  and  Lay,  1994;  Wu,  et  al. 
2002,  2003,  2004),  which  is  targeted  to  investigate  the  near-source  energy  partitioning  for 
regional  phases.  (2)  The  generalized  screen  propagator  (GSP)  method  (Wu,  1994,  1996;  Wu,  Jin 
and  Xie,  1996,  2000a,  b;  Xie  and  Wu,  2001;  Wu  and  Wu,  2001;  Wu,  et  al.,  1998,  1999,  2000). 
This  method  handles  laterally  varying  crustal  structures  including  volumetric  heterogeneities, 
intrinsic  attenuation  and  rough  free  surface.  Under  one-way  wave  equation  theory,  it  is  highly 
suited  for  investigating  high  frequency  and  long  distance  regional  wave  propagation.  (3)  The 
boundary  element  method  (Fu  and  Wu,  2000,  2001;  Fu,  et  al.,  2002;  Wu,  et  al.  2003).  This 
method  can  accurately  handle  uneven  surface  topography  and  interfaces  with  strong  impedance 
contrasts  and  is  suitable  for  investigating  interactions  between  waves  and  discontinuities.  These 
numerical  tools  teamed  together  enable  the  simulation  of  the  entire  regional  wave  excitation  and 
propagation  processes. 
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2.  LG-WAVE  EXCITATION  AND  NEAR-SOIJRCE  ENERGY  PARTITIONING 


2.1  Background 

For  explosion  sources,  there  are  major  questions  about  the  nature  of  excitation  of  S-wave 
dominated  phases  such  as  Lg,  and  there  are  similar  questions  regarding  the  relative  excitation 
effects  for  P/S-type  ratios  in  regional  phases,  particularly  given  the  huge  scatter  observed  in  both 
earthquake  and  explosion  data  populations.  Several  possible  near  source  energy  transfer 
mechanisms  have  been  proposed,  including  P-  Lg  scattering,  pS-Lg  conversion  at  the  free 
surface,  Rg-to-Lg  coupling,  S*-Lg  conversion,  spall  excitation  of  S,  tectonic  release  and  rock- 
damage  (e  g.,  Day  and  Mclaughlin,  1991;  Gupta  et  al.,  1992,  1997;  Wallace,  1991;  Gutowski,  et 
al.,  1984;  Vogfjord,  1997;  Johnson  and  Sammis,  2001).  There  is  continuing  controversy  about 
which  mechanism(s)  dominate  the  explosion  source  energy  partitioning  processes.  Recent 
observations  from  the  DOB  experiment  at  the  Shagan  Test  Site  provide  evidence  that  the  shallow 
source  exhibits  increased  S-wave  excitation,  which  supports  the  Rg-to-Lg  conversion  hypothesis 
(Myers,  et  al.,  2003;  Bonner,  et  al.,  2003).  However,  other  evidence  suggests  that  the  S-wave 
may  be  generated  primarily  by  the  nonspherical  part  of  the  explosion  source,  with  strong 
influence  from  the  free  surface  (Stevens,  et  al.,  2003).  Another  unsolved  puzzle  is  the  frequency 
dependence  observed  for  regional  P/S  discriminants.  Numerous  studies  (e  g.,  Walter  et  al.,  1995; 
Fisk  et  al.,  1996;  Taylor,  1996;  Taylor  and  Hartse,  1997;  Hartse  et  al.,  1997;  Fan  et  al.,  2002) 
have  demonstrated  that  regional  P/S  measurements  (Pn/Lg,  Pg/Lg,  Pn/Sn,  Pg/Sn)  provide  useful 
discrimination  of  earthquakes  and  explosions  at  frequencies  above  about  3  Hz,  but  the 
discrimination  performance  is  much  lower  at  low  frequencies.  There  is  very  incomplete 
theoretical  understanding  of  this  strong  frequency-dependent  phenomenon 

Early  numerical  simulations  investigated  various  propagation  properties  of  the  Lg  phase  using 
the  wavenumber  integral  method  (e  g.,  Bouchon,  1982;  Campillo  et  al.,  1984;  Campillo  and  Paul, 
1 992).  These  works  revealed  many  propagation  and  excitation  characteristics  of  regional  phases 
in  horizontally  stratified  crustal  models.  However,  with  a  1 D  model,  many  phenomena  linked  to 
the  laterally  varying  crustal  heterogeneity  cannot  be  generated  so  the  results  were  very  limited. 

Using  full  wave  numerical  techniques  such  as  finite-difference  (FD)  (Frankel  1986;  Vidale  and 
Helmberger,  1988;  Hayashi,  et  al.,  2001)  or  pseudo  spectrum  methods  (Kosloff  and  Baysal, 

1982;  Tessmer  and  Kosloff,  1994;  Orrey,  et  al.,  2003)  to  explore  seismic  wave  propagation  has 
the  advantage  that  complicated  two-  or  three-dimensional  structures  can  be  included  in  the  model 
and  the  complete  wave  field  can  be  calculated  These  numerical  methods  are  also  widely  used  for 
investigating  Lg-wave  excitation  and  propagation  in  complicated  crustal  waveguides.  Xie  and 
Lay  (1994)  investigated  Lg-wave  excitation  using  the  FD  method.  Jih  (1995,  1996)  investigated 
Rg-to-Lg  coupling  as  a  possible  Lg  excitation  mechanism.  Using  2D  and  3D  general  Fourier 
method,  Bonner,  et  al.  (2003)  investigated  Rg  and  Lg  generation,  and  partially  reproduced  the 
observed  spectrum  from  the  Depth  of  Burial  Experiment.  Stevens  et  al.  (2003)  investigated  the 
physical  basis  of  explosion  generated  S-waves  using  a  2D  nonlinear  FD  method,  which  handles 
axisymmetric  near  source  effects  including  spall,  cracking,  and  nonlinear  deformation. 

The  main  disadvantages  of  these  numerical  methods  that  provide  complete  synthetic 
seismograms  are  the  low  computation  efficiency  and  huge  computer  memory  requirement, 
especially  when  applied  to  investigate  the  characteristics  of  Lg  excitation.  For  the  purpose  of 
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small  nuclear  test  monitoring,  the  range  of  interest  for  Lg-wave  simulation  involves  a  broad 
frequency  band  (0.2  to  10  Hz)  and  long  propagation  distances  (up  to  1000  km  or  more).  At  the 
same  time,  the  factors  that  control  the  source  energy  partitioning  depend  on  the  detailed  source 
mechanism  and  fine  near  source  velocity  structure.  In  addition,  there  are  multiple  mechanisms 
that  may  potentially  contribute  to  the  energy  partitioning  process.  Numerous  parameters  need  to 
be  tested  to  investigate  the  characteristics  of  these  mechanisms,  especially  their  contributions  to 
the  frequency  dependent  features  of  observable  discriminants.  If  random  heterogeneities  are  to 
be  considered,  as  is  likely  to  be  important  for  high  frequency  signals,  the  results  have  to  be 
calculated  statistically  from  simulations  using  a  large  number  of  realizations.  This  limits  the 
approach  of  complete  FD  synthesis  for  actually  recording  geometries. 

Although  there  are  continuing  controversies  about  the  dominant  P-to-S  transfer  mechanisms 
affecting  regional  phases,  most  investigators  agree  that  appreciable  energy  from  explosion 
sources  is  converted  to  S- waves  in  the  near-source  region  (Myers,  2003).  The  physical  processes 
by  which  an  explosion  source  generates  regional  phases  can  be  described  as  energy  partitioning 
taking  place  in  the  near-source  region.  The  partitioned  energy  subsequently  propagates  through 
a  long  waveguide,  where  secondary  energy  partitioning  effects  may  take  place,  but  these  are  less 
associated  with  the  type  of  source  involved.  The  combined  processes  in  the  near-source  region 
and  along  the  path  naturally  separate  the  wave  field  energy  into  groups  of  observable  regional 
phases  according  to  their  slowness  and  group  velocity.  If  the  propagation  effect  is  not  the 
primary  goal  of  an  investigation,  there  is  a  strategy  to  avoid  calculating  the  immensely  time- 
consuming  long  distance  propagation  part  of  the  problem,  focusing  on  the  near-source  energy 
partitioning  effects.  Based  on  this  concept,  we  developed  a  method  based  on  the  finite-difference 
simulation  and  slowness  analysis  to  investigate  the  near-source  energy  partitioning  of  an 
explosion  source  This  method  investigates  the  partitioning  process  right  at  the  source  region. 

The  localized  analysis  thereby  provides  uncontaminated  information  isolating  the  physical 
processes  controlling  the  energy  partitioning. 

2.2  Methodology 

Many  authors  (e  g.  Frankel,  1989;  Xie  and  Lay,  1994;  Vogfjord,  1997)  have  pointed  out  that  for 
S-wave  energy  to  be  trapped  in  the  waveguide  reverberating  to  generate  the  Lg-wave,  it  must 
propagate  with  post  critical  angle  at  the  Moho  discontinuity.  However,  in  the  waveguide  and 
especially  in  the  near-source  region,  the  wavefield  is  highly  complicated.  It  is  impractical  to 
trace  each  phase  in  the  spatial-time  domain.  An  alternate,  but  equally  valid  way  of  tracking  the 
wave  energy  is  in  the  slowness  domain.  Multiply-reflected  waves  may  arrive  simultaneously  in 
time,  but  in  the  slowness  domain  their  energy  distribution  gives  clear  information  about  the  wave 
intensity,  slowness  and  propagation  direction  Several  methods  can  be  used  to  transfer  spatial¬ 
time  domain  data  into  slowness  (or  equivalently  wavenumber)  domain  information,  for  example, 
FK  analysis  or  beamforming  (slant  stacking).  Here  we  employ  the  slant  stack  method,  but  with 
the  novelty  that  we  use  localized  slowness  analysis,  working  simultaneously  on  both  space  and 
slowness  domain.  Based  on  this  method  (introduced  by  Xie  and  Lay,  1 994  and  further 
developed  by  Wu  et  al ,  2002,  2003),  the  wave  field  at  short  distance  is  sampled  with  different 
spatial,  time  and  frequency  windows.  Then  the  energy  is  characterized  as  a  function  of  phase 
velocity  based  on  the  slowness  analysis.  The  output  provides  localized  information  in  time. 
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space,  frequency  and  slowness  domains,  which  provides  the  overall  information  for  investigating 
the  characteristics  of  energy  partitioning. 


Figure  2. 1  illustrates  the  basic  concept  of  the  finite-difference  simulation  plus  slowness  analysis 
(FDSA)  method  for  investigating  the  near-source  Lg-wave  excitation  and  energy  partitioning.  In 
Figure  2.1,  (A)  shows  the  finite-difference  simulation  snap-shot  extended  out  to  regional 
distances  where  internal  waveguide  slowness  analysis  can  be  conducted  as  well  as  where  surface 
theoretical  seismograms  can  be  computed;  (B)  is  the  near  source  model  composed  of  different 
types  of  source  representations  and  fine  near  source  structures  that  can  be  used  to  test  different 
energy  partitioning  mechanisms;  (C)  is  a  summary  of  the  slowness  analysis  which  sorts  the 
energy  from  the  short  distance  array  data  and  predicts  the  trapped  waveguide  energy  at  larger 
distances. 


fine  near  source 
structure 


short  distance  array  observation 
inside  the  waveguide 


source 


Moho 


distance 
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A.  finite-difference  simulation 


Figure  2.1.  Cartoon  showing  the  finite-difference  simulation  and  slowness  analysis  method.  (A)  The  finite- 
difference  modeling  used  to  propagate  the  wavefield  w  ithin  the  crustal  waveguide  is  illustrated  by  a 
snap-shot  of  the  displacement  field.  (B)  Fine  near  source  velocity  model  used  to  investigate  near 
source  energy  partitioning.  (C)  Energy  flux  at  the  short  distance  is  calculated  and  used  to  predict  the 
long  distance  observations.  These  predictions  arc  corroborated  by  analysis  of  surface  sythctics  at 
large  distances. 

Figure  2.2  shows  different  energy  distributions  for  slowness  analysis.  The  energy  can  either  be 
expressed  as  a  function  in  slowness  domain,  or  as  a  function  in  mixed  slowness-depth  domain. 
The  measurement  can  be  made  in  the  entire  cross  section  of  the  waveguide.  The  upper  mantle  S- 
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wave  slowness  \lvs.mantie  is  marked  for  different  cases.  All  S-wave  energy  (either  directly 
radiated  from  the  source  or  converted  from  other  wave  types)  with  apparent  horizontal  slowness 
larger  than  the  upper  mantle  S-wave  slowness  will  be  trapped  in  the  crustal  waveguide  and  will 
contribute  to  the  guided  S  phase  unless  subsequent  scattering  causes  it  to  leak  out.  Effectively, 
the  relative  amplitude  of  the  trapped  S-wave  energy  can  be  determined  at  a  short  distance  The 
strength  of  these  S-waves  will  depend  on  two  factors,  the  intensity  of  the  S  energy  radiated  or 
converted  in  the  near-source  region,  as  well  as  the  percentage  of  the  energy  that  can  intrinsically 
be  trapped  in  the  waveguide.  The  slowness  analysis  is  done  within  spatial,  time  and  frequency 
windows.  The  localized  slowness  analysis  allows  the  energy  flux  in  the  waveguide  to  be 
investigated  in  mixed  domains  of  space,  travel-time,  slowness,  and  frequency.  The  energy 
distribution  in  these  joint  domains  has  a  better  chance  of  being  separated  than  in  conventional 
seismogram  analyses,  and  the  underlying  physical  mechanisms  controlling  the  energy 
partitioning  can  thus  be  investigated. 
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Figure  2.2.  Slowness  domain  display.  Left:  energy  distribution  in  2D  slowness  domain,  right:  energy 

distribution  in  mixed  horizontal  slowness-depth  domain.  Vertical  lines  indicate  the  upper-mantle  S- 
wave  slowness.  Energy  falling  on  the  right-hand  side  of  these  lines  will  he  trapped  in  the  waveguide 
and  will  form  guided  S  phases  like  Lg. 


2.3  Verification  of  the  method 

To  verify  our  method,  we  computed  a  suite  of  large  2D  FD  simulations  (see  Figure  2. 1 ).  For  all 
numerical  examples  calculated  in  this  report,  we  will  use  the  EK  (Eastern  Kazakh)  model 
(Priestley  et  al.,  1988)  as  the  background  and  modify  it  by  adding  different  laterally  varying 
structures.  A  range  of  sources  and  near-source  structures  were  used  in  the  simulations.  Surface 
receivers  at  long  distance  are  computed  to  provide  conventional  regional  phase  synthetics  that 
we  treat  as  distant  observations  of  the  wave  field.  Vertical  arrays  embedded  in  the  crustal 
waveguide  at  short  distances  are  used  for  the  slowness  analysis. 


Figure  2.3  shows  the  results  for  a  cross  section  at  180  km  in  the  crustal  waveguide.  In  each 
panel,  the  vertical  axis  is  the  horizontal  slowness.  The  horizontal  axis  is  time.  The  group  velocity 
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is  also  marked  on  the  top  of  each  panel.  The  dashed  line  marks  the  upper  mantle  S-wave 
velocity.  The  solid  circles  are  energy  picked  from  the  slowness  domain  and  summed  up  for  the 
entire  cross  section.  The  left  column  is  for  low  frequency  (0.3  -1.5  Hz)  and  the  right  column  is 
for  high-frequency  (2. 0-5.0  Hz).  The  top  panel  is  for  background  velocity  model.  It  is  important 
to  note  that  the  Rg  mainly  appears  as  a  low-frequency  signal  and  the  Pg  mainly  appears  as  the 
high-frequency  signal.  We  first  focus  on  the  low  frequency  results.  From  the  top  panel  we  can 
see  Pg-wave  and  strong  fundamental  mode  Rg-wave.  However,  such  a  high  velocity  crust 
generates  very  poor  Lg- waves.  Within  the  Lg  wave  window,  only  very  weak  energy  can  be 
observed  which  may  result  from  the  S*-wave.  In  the  second  panel,  a  random  velocity  patch 
located  at  depth  0.0-2. 5  km  is  added  to  the  background  model  (refer  to  Figure  2.1).  Compared 
with  the  background  model,  through  scattering,  the  Rg-wave  loses  considerable  energy  while  the 
Lg-wave  gains  energy.  Panels  3  to  5  are  similar  to  panel  2,  except  random  patches  are  located  at 
depths  2. 5-5.0  km,  5. 0-7.5  km  and  10.0-12.5  km,  respectively.  The  tendency  is:  the  shallower  the 
random  patch  the  more  energy  is  transferred.  This  also  provides  additional  evidence  that  the  low 
frequency  Lg  energy  may  come  from  the  Rg-wave,  since  the  Rg  energy  is  mostly  concentrated  at 
the  shallow  depth 

Figure  2.4  shows  the  comparison  between  waveguide  energy  flux  at  180  km  calculated  using 
slowness  analysis  (left  panel)  and  at  450  km  calculated  from  the  surface  measurement  using 
conventional  Lg  wave  processing  technique  (right  panel).  For  simplicity,  we  label  the 
predictions  of  trapped  energy  as  “Lg”  energy.  Several  random  velocity  patches  are  incorporated 
in  the  near-source  region  to  provide  the  wave  field  coupling.  The  vertical  extent  of  the  patches  is 
2.5  km,  with  RMS  velocity  fluctuation  of  10%  and  the  positioning  of  the  patches  varies  from  0  to 
1 5  km  in  depth.  The  results  show  the  relative  energy  changes  (short  color  bars)  of  Pg-,  Lg-  and 
Rg-waves,  as  functions  of  the  depth  range  of  the  random  velocity  scattering  patch  for  each  case. 
The  thin  lines  indicate  the  reference  energy  level  for  a  velocity  model  without  any  random 
velocity  scattering  patch.  Short  color  bars  indicate  the  energy  changes  due  to  the  near  source 
scattering  at  different  depths.  These  calculations  demonstrate  that  the  relative  changes  of  energy 
for  each  of  these  phases  obtained  in  the  waveguide  slowness  analysis  correspond  very  closely  to 
those  obtained  on  the  free  surface  at  long  distance,  even  for  dramatic  changes  such  as  the  strong 
scattering  of  Rg  which  occurs  when  the  scattering  patch  is  just  below  the  surface.  These 
calculations  confirm  that  the  near-source  slowness  analysis  within  the  waveguide  correctly 
predicts  the  surface  regional  observations  at  greater  distance;  this  is  the  heart  of  our  strategy. 

To  examine  the  energy  flux  measurements  obtained  at  even  shorter  distances,  we  compared  the 
slowness  analysis  measurements  at  1 80  km  with  measurements  at  60  km  distance.  Figure  2.5 
gives  examples  for  varying  phase  velocities  at  these  two  distances.  The  wave  fields  at  these 
distances  show  quite  different  features  in  the  slowness-depth  domain  but  give  consistent 
predictions  for  the  wave  field  energy  fluxes.  Figure  2.6  compares  the  corresponding  wave  guide 
energy  measured  at  60  km  and  1 80  km  distances.  Although  a  wide  range  of  near-source 
structures  and  source  depths  are  used  to  generate  these  measurements,  they  show  very  good 
linear  relationships.  This  further  verifies  that  we  can  use  a  small  model  to  predict  the  energy 
distribution  at  greater  distances. 
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Figure  2.3.  Slowness  analysis  for  velocity  models  with  random  patches  located  at  different  depths.  The  left 
column  is  for  the  low-frequency  hand  (0.3-1. 5  Hz)  and  the  right  column  is  for  the  high-frequency 
band  (2.0-5.0  Hz).  The  top  panels  show  the  background  velocity  model,  and  the  rest  of  panels  show 
random  patches  at  different  depths.  For  details,  see  the  text. 
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Figure  2.4.  Comparison  between  waveguide  energy  flux  at  180  km  (left  panel)  and  the  wave  energy  on  the 

surface  at  450  km  (right  panel)  for  Pg,  Lg  and  Rg  windows.  Different  near  source  velocity’  models  arc- 
used  in  the  calculation.  The  frequency  range  is  0.3-1. 2  Hz.  Shown  in  each  small  panel  is  relative 
energy,  the  horizontal  axis  is  the  depth  of  the  random  patches,  the  vertical  axis  is  the  relative  energy  , 
thin  line  indicates  the  energy  level  for  background  model,  short  color  bars  indicate  the  energy 
changes  due  to  the  near  source  scatterings. 
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Figure  2.5.  Examples  showing  the  slowness  analyses  at  distances  of  60  km  (upper  panel)  and  180  km  (lower 
panel).  The  energy  that  can  be  trapped  is  indicated  in  the  figure.  The  wavefields  show  quite  different 
features  at  these  distances  but  give  consistent  predictions  for  the  trapped  energy  fluxes. 


Energy  flux  at  60  km  Energy  flux  at  60  km  Energy  flux  at  60  km 


Figure  2.6.  Comparison  between  trapped  waveguide  energy  measured  at  60  km  (horizontal  coordinate)  and 
180  km  (vertical  coordinate)  for  frequency  hands  0.3-1.2  Hz,  1.0-2.5  Hz  and  2.0-5.0  Hz.  The  results 
are  from  different  velocity  models  and  source  depths.  Different  colors  represent  different  velocity 
models.  The  results  show  a  general  linear  relationship  for  all  frequency  hands. 
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Figure*  2.7.  Summary  of  the*  P-pS-Lg  wave  coupling.  From  left  to  right  are  velocity  models,  wavefield 
snapshots  and  energy  distributions  in  the  slowness-depth  domain.  Details  sec  the  text. 


2.4.  Investigating  the  Regional  Phase  Excitation 


2.4. 1  The  P-Lg  am I pS-Lg  conversion 

We  first  assess  the  importance  of  near-source  lateral  velocity  variations  that  cause  coupling 
between  P,  pS,  and  Lg.  Figure  2.7  summarizes  how  the  depth  of  heterogeneity  affects  the  P-to-S 
coupling.  From  left  to  right  are  velocity  models,  wave  field  snapshots  and  energy  distributions 
in  the  slowness  domain.  Slowness  analysis  is  used  at  a  distance  of  16  km  and  depth  0-10  km 
We  use  different  source/model  combinations  in  the  simulations  and  test  their  effects  on  the  P-to- 
S  conversion.  This  example  shows  that  without  the  lateral  velocity  heterogeneity,  the  energy  of 
the  free-surface  reflected  pS  wave  falls  to  the  left  of  the  upper  mantle  S-wave  slowness  and 
therefore  cannot  form  trapped  energy  at  regional  distance.  With  a  random  velocity  layer  added  to 
the  model,  the  pS-wave  front  becomes  distorted  and  its  slowness  content  crosses  the  upper 
mantle  S-slowness.  Assessing  the  energy  that  builds  to  the  right  of  the  upper  mantle  S-slowness 
allows  us  to  predict  the  Lg-wave  component  associated  with  P-to-S  conversion.  As  can  be  seen 
from  numerical  examples,  the  depth  of  the  source  and  the  depth  of  the  scattering  region  have 
strong  effects  on  the  P-to-S  conversion  and  the  partitioning  of  energy  into  trapped  and  leaking 
signals. 


2.4.2  The  Rg-to-Lg  coupling 

Rg-to-Lg  coupling  has  long  been  proposed  as  a  potential  mechanism  for  low  frequency  Lg-wave 
excitation.  However,  there  is  disagreement  as  to  whether  this  is  a  dominant  mechanism  or  not. 
This  is  mainly  due  to  the  lack  of  knowledge  about  the  actual  mechanism  of  this  process,  for 
example,  where  this  process  happens,  how  much  energy  can  be  transferred  from  Rg  to  Lg,  where 
the  balance  of  the  energy  goes,  and  what  is  the  frequency  dependence  of  this  process?  Although 
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it  is  generally  agreed  that  Rg-to-Lg  coupling  is  caused  by  scattering  in  the  near-surface 
environment,  systematic  investigation  of  this  process  is  still  lacking.  Since  the  coupling  process 
links  the  high  end  of  the  Rg  band  to  the  low  end  of  the  Lg  band,  the  statistical  feature  of  the 
lateral  velocity  or  topography  variation  must  play  an  important  role  in  selecting  the  coupling. 
Any  numerical  test  targeted  to  solve  this  problem  must  comprehensively  take  care  of  the 
excitation  of  the  Rg  (the  source  depth,  shallow  structure,  CLVD  component  of  the  source),  the 
frequency  dependent  scattering  process  (shallow  structures  or  rough  free  surface  with  different 
statistical  parameters)  and  sorting  of  the  scattered  energy  to  find  the  component  eventually 
trapped  in  the  waveguide  as  Lg  energy.  Using  the  FDSA  method,  we  have  conducted  some 
numerical  investigations  of  this  problem. 
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Figure  2.8.  The  apparent  relationship  between  the  Lg-  and  Rg-wave  energy  at  180  km  in  the  waveguide  vs. 

the  depth  of  near  source  random  patches.  Frequency  range  is  0.3-1. 2  Hz.  Red  dash  lines  are  for  the 
background  layered  model.  Blue  bars  are  for  models  with  near  source  random  patches  at  different 
depths. 


We  conducted  two  types  of  tests:  a)  for  a  shallow  explosion  source,  we  locate  scattering  regions 
at  different  depths;  and  b)  for  a  shallow  scattering  region,  sources  are  located  at  different  depths. 
Figure  2.8  shows  the  waveguide  energy  fluxes  for  Lg  and  Rg  versus  the  depth  of  near-source 
heterogeneities.  The  frequency  range  is  0.3- 1.2  Hz,  the  overlapping  domain  for  the  two  phases. 
Red  lines  are  for  the  background  model.  Blue  bars  are  for  models  with  near-source  random 
layers  at  different  depths.  For  models  with  shallow  scattering,  Rg  loses  energy  while  Lg  gains 
energy.  Figure  2.9  gives  the  apparent  relationship  between  the  source  depth  and  the  Rg  to  Lg 
coupling.  The  three  frequency  bands  are  0.3-1. 5,  1. 0-2.0  and  2.0-5. 0  Hz,  respectively.  The 
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sources  are  located  at  different  depths  and  scattering  is  provided  by  adding  a  near  source  random 
patch  at  the  top  of  the  crust.  A  shallower  source  raises  both  Rg  and  Lg  energy  while  shallow 
scattering  increases  the  Lg  energy  but  decreases  the  Rg  energy.  This  is  consistent  with  some 
observations.  The  apparent  Rg-Lg  coupling  occurs  in  the  low  frequency  band  and  Pg-Lg 
coupling  occurs  in  the  high  frequency  band.  Although  these  preliminary  results  provide  an 
apparent  relationship  similar  to  some  observations,  they  do  not  yet  provide  robust  evidence  that 
the  energy  lost  from  Rg  is  necessarily  transferred  to  the  Lg,  since  the  S*-wave  also  has  a  similar 
excitation  function,  i.e.,  it  favors  shallow  source  and  lower  frequency  excitation 
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Figure  2.9.  Effects  of  shallow  near  source  scattering  on  the  energy  partitioning  of  Pg-,  Lg-  and  Rg- waves. 

Sources  are  located  at  different  depths  between  0.5  and  10  km.  Results  for  the  background  model  are 
shown  with  red  bars  and  results  for  models  with  random  patches  are  shown  with  blue  bars. 


2.4.3  The  S*-to-Lg  excitation  spectrum 

For  shallow  spherical  explosions,  depending  on  the  situation,  the  S*-wave  may  become  a 
significant  contributor  to  Lg  (Gutowski,  et  al.,  1984;  Xie  and  Lay,  1994;  Vogfjord,  1997).  The 
amplitude  of  S*  can  be  large  as  long  as  the  source  depth  is  within  a  fraction  of  a  wave  length 
from  an  interface.  This  makes  its  excitation  highly  dependent  on  the  source  depth  and  frequency. 
Figure  2. 10  shows  snapshots  for  explosion  sources  at  0.5  km  and  3.0  km,  respectively.  This 
clearly  shows  that  a  shallow  source  generates  larger  S*-  and  Rg-waves.  With  the  FDSA  method, 
their  energy  can  be  easily  investigated  near  the  source  region  because  of  the  distinct  slowness  for 
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these  phases.  Figure  2. 1 1  gives  slowness  analysis  at  a  60  km  distance  for  explosions  at  different 
depths.  We  can  isolate  and  quantify  the  S*  energy  even  within  the  complicated  wavefield,  which 
is  very  difficult  to  do  with  isolated  surface  synthetics.  Using  a  band  pass  filter,  the  excitation 
function  of  the  S*-wave  can  be  obtained.  This  can  be  compared  with  various  observations  such 
as  the  spectral  content  of  Sn  and  Lg,  and  will  allow  us  to  assess  whether  S*  is  in  fact  an  effective 
Lg  source  relative  to  other  mechanisms. 
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Figure  2.10.  Wavefield  snapshots  for  explosion  sources  at  different  depths.  The  source  depths  for  left  and 

right  panels  are  0.5  km  and  3.0  km,  respectively.  A  shallower  explosion  is  a  more  efficient  source  for 
S*  and  Rg- waves. 


We  have  shown  some  results  that  have  intrinsic  frequency  dependence.  These  frequency 
variations  are  rooted  in  the  underlying  physical  processes  and  are  usually  controlled  by  different 
characteristic  scales.  For  example,  the  excitation  of  Lg  by  S’",  Rg-to-S  scattering  and  spall  are 
all  highly  source  depth  dependent.  The  excitation  spectra  of  specific  phases  depict  the  frequency 
dependence  of  these  processes.  Figure  2. 12  gives  the  S*-to-Lg  excitation  spectra  for  sources  in 
different  models  and  at  different  depths  The  result  clearly  shows  that  the  excitation  is  generally 
favorable  to  lower  frequency  and  shallow  source  depth.  However,  it  is  also  model  dependent. 
For  a  model  with  homogeneous  upper  crust  (left  panel),  the  distribution  has  a  simple 
monotonous  tendency.  While  for  the  EK  model,  which  has  an  interface  at  depth  1 .0  km  (right 
panel),  the  excitation  spectrum  has  a  maximum  at  depth  1  km  and  a  more  complicated  frequency 
dependency.  This  may  result  from  the  particular  velocity  model  used.  The  frequency 
dependence  of  these  excitation  spectra  can  provide  the  basis  for  evaluating  the  dominant 
mechanisms  for  Lg-wave  excitation,  or  for  down-weighting  particular  mechanisms  as  being  less 
significant.  It  will  also  provide  the  relationship  between  the  observations  and  the  characteristics 
of  the  source  and  near  source  structure 
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Figure  2.11.  Slowness  analysis  for  investigating  the  S*-wave  energy  eontributing  to  the  Lg-wave.  Different 
rows  are  for  different  source  depths.  Dash  line  rectangles  indicate  the  time-space- slow  ness  windows 
used  to  pick  the  S*  energy  . 


Figure  2.1 1.  S*-to-Lg  excitation  spectra  for  sources  at  different  depths.  Left:  a  model  with  homogeneous 
upper  crust,  and  Right:  EK  model.  Details  see  the  text. 
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2.4.4  The  frequency  dependent  discriminants 

Frequency  dependent  P/S  ratios  will  depend  on  the  excitation  functions  of  multiple  phases. 

Figure  2.12  shows  observed  Log  (Pn/Lg)  amplitude  ratios  in  the  frequency  band  between  0.5  and 
10  0  Hz  for  explosions  and  earthquakes  (Fisk  2003).  The  two  populations  exhibit  poor 
separation  below  about  2  to  3  Hz  and  good  separation  at  higher  frequencies.  Figure  2. 1 3  gives  a 
simulation  of  the  Log(P/Lg)  amplitude  ratio  in  the  frequency  band  0.2  to  5.2  Hz.  Two  types  of 
sources  are  used,  an  isotropic  explosive  source  at  0.5  km  depth  representing  a  typical 
“explosion”,  and  a  45  degree  dip-slip  dislocation  source  at  3.0  km  depth  representing  a  typical 
“shallow  earthquake”.  Two  velocity  models  are  used  in  calculation,  EK  and  EK  plus  a  shallow 
near  source  random  layer  (depth  0-2.5  km,  distance  5-25  km).  The  slowness  analysis  is 
conducted  at  a  distance  of  60  km.  For  the  explosive  source  in  the  EK  model,  the  Lg  energy 
mainly  comes  from  the  S*  phase.  For  the  EK  plus  random  layer  model,  it  is  expected  that  S*,  P- 
to-Lg  and  Rg-to-Lg  couplings  will  jointly  contribute  to  the  Lg  energy.  Compared  with  the 
explosive  source,  the  somewhat  deeper  dislocation  source  is  less  affected  by  the  shallow  random 
layer.  The  result  shows  a  general  feature,  which  is  similar  to  that  from  the  observations  (Figure 
2. 12).  At  the  low  frequency  end,  the  earthquakes  and  explosions  have  similar  P/Lg  ratios  but  at 
higher  frequencies,  the  explosive  source  generates  less  S  energy. 


Figure  2.12.  Log  Pn/Lg  amplitude  ratios  in  frequency  bands  from  0.5  to  10.0  Hz  for  20  Balapan  nuclear 
explosions  (asterisks)  and  24  earthquakes  (squares)  recorded  by  WMQ  (Fisk,  2003).  The  two 
populations  exhibit  poor  separation  below  about  2  to  3  Hz.  and  good  separation  at  higher  frequences. 
SNRs  for  three  of  the  nuclear  explosions  are  marginal  at  frequencies  above  ~8  Hz 
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It  is  expected  that  with  more  complexities  being  added  to  the  model  (CLVD,  different 
overburden  structures  and  scattering  models,  source  depth  changes,  different  source  time 
functions,  etc  ),  these  frequency  dependent  features  will  be  complicated,  but  they  can  be 
systematically  mapped  out.  This  is  the  essential  undertaking  for  placing  regional  phase  energy 
partitioning  on  a  sound  physical  basis.  The  numerical  modeling  will  set  up  a  link  between  the 
physical  model  and  observable  frequency  dependent  features,  which  can  provide  guidelines  for 
improving  the  empirical  discrimination  relations.  The  feedback  from  the  observations  will  guide 
the  simulations  towards  the  correct  physical  models  as  well. 
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Figure  2.13.  Simulated  Log  (P/Lg)  amplitude  ratios  in  frequency  domain.  The  circle,  diamond,  triangle  and 
cross  represent  explosion  source  in  the  EK  model,  explosion  source  in  the  EK  model  with  shallow 
random  patch,  dislocation  source  in  the  EK  model  and  dislocation  source  in  the  EK  model  with 
shallow  random  patch.  The  explosion  source  is  located  at  depth  0.5  km  and  the  dislocation  source  is 
located  at  depth  3.0  km. 
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3.  THE  PATH  EFFECT  ON  LG  WAVE  PROPAGATION 


3.1  Background 

High-frequency  regional  wave  propagation  in  complex  crustal  waveguides  is  one  of  the  most 
challenging  problems  in  theoretical  and  computational  seismology.  For  the  applications  of 
regional  waves  to  various  geophysical  problems,  a  good  understanding  of  their  propagation, 
scattering,  attenuation  and  wave-type  conversion,  as  well  as  the  availability  of  tools  to  simulate 
and  analyze  these  phenomena  in  complex  crustal  structures  (including  rough  topography,  Moho 
undulation  and  volumetric  velocity  heterogeneities  of  different  scales)  are  crucial  Regional  wave 
tomography  for  crustal  structures,  path  correction  for  discrimination  and  yield  estimation  of  low- 
yield  nuclear  tests,  location  determination  of  earthquakes  or  underground  explosions  using 
regional  phases  are  examples  of  the  possible  applications.  Monitoring  the  CTBT 
(Comprehensive  Nuclear-Test-Ban  Treaty)  at  regional  distances  is  even  more  demanding  of  the 
simulation  and  analyzing  tools.  For  this  purpose,  simulation  algorithms  are  desirable  for 
generating  synthetic  waveforms  for  high  frequencies  and  propagation  distances  greater  than  1 000 
km. 

Substantial  efforts  have  been  made  in  modeling  regional  wave  propagation.  Methods  based  on 
layered  earth  models,  such  as  the  reflectivity  and  mode  summation  methods  (e  g.,  Bouchon  et  al., 
1985;  Kennett,  1989,  1990;Maupin,  1989;  Baumgardt,  1990;  Campillo,  1990;  Campillo  and 
Paul,  1992;  Campillo  et  al.,  1993;  Gibson  et  al.,  1994)  have  very  high  efficiency  and  can  be 
applied  to  relatively  high  frequencies,  but  they  can  be  used  only  for  very  simplified  cases  with 
layered  or  smoothly  varying  layered  models.  Modeling  techniques  that  can  treat  realistic  3D 
heterogeneous  media  rather  than  smoothly  varying  layered  media  are  needed  to  test  and  study 
many  observations  and  hypotheses.  Sudden  changes  of  crustal  thickness,  strong  lateral  variations 
and  irregular  3D  heterogeneities  are  among  the  problems  requiring  new  modeling  methods.  As 
pointed  out  by  Campillo  et  al.  (1993),  actual  Lg  amplitudes  were  reduced  more  than  10  times  for 
paths  passing  through  an  anomalous  zone  on  the  east  side  of  the  Alpine  range,  while  the 
modeling  results  using  existing  methods  (including  the  effect  of  known  large-scale  lateral 
structural  variation)  only  account  for  20  -  30  %  of  the  amplitude  reduction.  Other  mechanisms 
such  as  scattering  and  attenuation  by  small-scale  heterogeneities  must  be  taken  into  account. 

Kennett  ( 1 984,  1 998),  Maupin  and  Kennett  ( 1 987)  developed  a  coupled  mode  method  for 
calculating  guided  seismic  waves  in  horizontally  varying  structures.  The  method  works  well  for 
relatively  low  frequency  waves  in  moderately  heterogeneous  models.  However,  the 
implementation  of  the  method  for  high  frequency  3-D  models  still  requires  formidable 
computational  efforts. 

Cormier  and  Anderson  (1996,  1997)  applied  elastic  Bom  scattering  (in  the  regime  of  Rayleigh 
scattering)  to  the  locked-mode  solution  for  plane  layered  media  to  calculate  the  effects  of  small- 
scale  heterogeneities.  However,  the  approximation  is  limited  to  single  scattering  and  is  only  good 
for  heterogeneities  with  scales  much  smaller  than  the  wavelength.  Ray  methods  have  very 
limited  success  in  modeling  regional  waves  due  to  the  chaotic  behavior  of  rays  caused  by  the 
multiple  reflections  of  the  ffee-surface  and  Moho.  Keers  et  al.  (1996a,  b)  applied  the  Maslov 
integral  method  to  avoid  the  caustics  and  pseudo-caustics  (caustics  of  plane  waves)  by  working 
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in  the  phase-space.  However,  when  chaos  develops  in  the  ray  system,  there  are  more  complicated 
caustics  for  which  the  Maslov  method  does  not  work.  In  addition,  ray-tracing  computation  is 
very  time-consuming  in  this  case.  An  alternative  and  flexible  approach  using  ray-tracing  has 
been  developed  by  Kennett  (1986),  Bostock  and  Kennett  (1990)  and  Kennett  et  al.  (1990),  in 
which  ray  diagrams  are  used  to  study  Lg  waves  when  crossing  structural  boundaries.  The  method 
agrees  well  with  modal  calculations  and  can  be  applied  to  surface  topography,  3-D  crustal 
structures  and  other  cases.  However,  the  method  cannot  provide  information  on  wave 
phenomena  for  complicated  waveguides. 

Finite-difference  methods  are  general  numerical  methods  and  have  been  extensively  applied  to 
regional  wave  propagations  (e  g.,  Xie  and  Lay,  1994;  App  et  al.,  1996;  Goldstein  et  al.,  1996, 
1997,  1999;  Husebye  and  Ruud,  1996;  Jih,  1996;  Nolte  et  al.,  1996;  Jones  et  al.,  1997; 
McLaughlin  and  Wilkins,  1997;  Bradley  and  Jones,  1998,  1999)  and  pseudo-spectral  methods 
(e  g.,  Kosloff  et  al.,  1990;  Orrey  et  al.,  2003;  Archambeau  et  al.,  1996;  Schatzman,  1996; 
Furumura  and  Kennett,  1997).  Theoretically  the  method  can  deal  with  arbitrarily  heterogeneous 
media.  However,  it  is  necessary  to  use  very  dense  spatial  sampling  to  avoid  grid  dispersion  for 
long  distance  regional  wave  propagation.  The  capability  of  the  present  day  computers  usually 
restricts  them  to  short  propagation  ranges  and  relatively  low  frequencies,  which  prevents  them 
from  being  applied  to  more  realistic  cases. 

The  state  of  the  art  of  the  traditional  simulation  techniques  for  regional  waves  has  its  application 
to  relatively  low  frequencies  and  short  propagation  distances.  Correspondingly,  the  volume 
heterogeneities  and  surface  irregularity  in  the  crustal  models  are  limited  to  rather  large  scales. 
However,  high  frequency  regional  waves  up  to  20  Hz  and  higher  have  been  observed  over 
distances  ranging  from  a  few  hundred  kilometers  to  more  than  one  thousand  kilometers  (e  g.,  Ni 
et  al,  1996;  Herman  et  al.,  1997;  Lay  et  al.,  1999).  Since  high-frequency  waves  can  be  used  for 
event  locations  with  high  accuracy,  simulation  and  modeling  of  high-frequency  regional  wave 
propagation  are  very  desirable  for  many  applications.  For  high  frequency  wave  propagation, 
scattering  and  attenuation,  the  role  of  small-scale  heterogeneities  and  surface  roughness  is 
important. 

The  existence  of  small-scale  heterogeneities  in  the  crust  and  the  associated  seismic  wave 
scattering  has  long  been  known  among  seismologists  (e  g.,  Aki,  1980;  Wu  and  Aki,  1988,  1989, 
1990;  Sato  and  Fehler,  1998).  However,  the  effects  of  these  heterogeneities  on  guided  wave  (Lg) 
propagation  in  the  crust  have  not  been  seriously  explored.  The  reasons  may  be  the  following. 
First,  the  spectra,  strength  and  distribution  of  the  small-scale  heterogeneities  in  different  regions 
are  not  clear  Very  few  data  sets  can  be  used  to  characterize  the  paths  concerned.  Second,  there  is 
a  lack  of  analytical  and  numerical  tools  to  model  or  analyze  their  influence  on  the  guided  wave 
propagation.  The  theory  of  wave  propagation  in  unbounded  random  media  has  been  well 
developed  However,  for  waves  in  complex  crustal  waveguides  with  random  heterogeneities,  the 
theoretical  difficulties  are  overwhelming,  and  no  analytical  tool  is  available  for  performing 
realistic  calculations.  Therefore,  numerical  methods  for  simulating  regional  wave  propagation  in 
complex  waveguides  with  small-scale  heterogeneities  are  highly  desirable.  It  has  become  clear 
that  small-scale  heterogeneities  are  widely  distributed  in  tectonically  active  regions.  Strong 
topographic  variation  is  the  manifestation  of  tectonically  active  regions  and  often  the  indication 
of  small-scale  heterogeneities.  Figure  3. 1  gives  a  topographic  profile  (top  panel)  and  its  power 
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spectrum  (bottom  panel)  for  a  path  crossing  Tibet  Plateau  The  slope  of  the  spectral  roll-off  is 
close  to  1/A,  a  flicker  noise  spectrum,  very  rich  in  small-scale  variations.  This  spectrum  is  similar 
to  the  observations  of  the  sonic  well-log  in  the  KTB  super-deep  continental  drilling  well  (Wu  et 
al.,  1994;  Jones  and  Holliger,  1997;  Goff  and  Holliger,  2002).  The  spectra  there  have  also  a  1/A 
slope.  Recently,  Goff  and  Holliger  (2002)  explained  the  1/A  spectra  as  a  combination  of 
hierarchical,  multi-scale  heterogeneities.  Overall,  the  1/A  spectra  demonstrate  richness  of  small- 
scale  heterogeneities. 
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Figure  3.1  The  topographic  profile  and  its  power  spectrum  for  a  path  crossing  the  Tibetan  Plateau. 

Recently,  the  generalized  screen  method  has  been  introduced  into  seismic  wave  simulations  and 
applied  to  the  problems  of  both  exploration  and  theoretical  seismology.  The  generalized  screen 
method  is  based  on  the  one-way  wave  equation  and  the  one-return  approximation.  The  one-way 
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generalized  screen  propagator  (GSP)  neglects  backscattered  waves,  but  correctly  handles  all  the 
forward  multiple-scattering  effects,  e  g.,  focusing/defocusing,  diffraction,  interference,  and 
conversion  between  different  wave  types.  The  one-retum  approximation  is  also  called  the  De 
Wolf  approximation  (De  Wolf,  1971,  1985),  which  neglects  the  reverberation  between  screens 
and  can  simulate  multiple-forescattenng-single-backscattering  (MFSB).  Significant  progress  has 
been  made  on  the  development  of  an  elastic  complex  screen  (ECS)  method  for  modeling  elastic 
wave  propagation  and  scattering  in  arbitrarily  complicated  structures  (Wu,  1994,  1996;  Xie  and 
Wu,  1995;  Wild  and  Hudson,  1998;  Wu  and  Wu,  2001).  The  method  is  two  to  three  orders  of 
magnitude  faster  than  the  elastic  finite-difference  method  for  medium  size  3D  problems.  The 
screen  method  has  been  successfully  used  in  forward  modeling  (Wu,  1994;  Wu  and  Huang, 

1995;  Xie  and  Wu  1995,  1996,  1999,  2001;  Wu  and  Wu,  2001)  and  as  backpropagators  for 
seismic  wave  imaging/migration  in  both  acoustic  and  elastic  media  (e  g.,  Wu  and  Xie  1994; 
Huang  and  Wu  1996;  Xie  and  Wu,  1998,  2004;  Huang  et  al.,  1999a,b;  Jin  and  Wu,  1999;  Jin  et 
al„  1999). 
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Figure  3.2  (A)  Upper  panel:  geometry  using  the  screen  method  to  simulate  Lg  wave,  and  (B)  lower 
panel:  sketch  show  ing  the  interaction  between  the  incident  waves  and  the  thin  slab. 
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3.2  The  Generalized  Screen  Propagator  for  Guided  Waves 

In  the  crustal  waveguide  environment,  major  wave  energy  is  carried  by  forward 
propagating  waves  bouncing  up  and  down  between  the  free  surface  and  major  geophysical 
discontinuities  such  as  the  Moho  and  Conrad  discontinuities.  Beyond  the  critical  angle,  these 
waves  are  systematically  dominated  by  small  angle  waves  trapped  in  the  crustal  waveguides. 
Therefore,  the  neglect  of  backscattered  waves  in  the  propagation  will  not  distort  the  main 
features  of  regional  waves  in  most  cases.  With  this  approximation,  the  method  becomes  a 
forward  marching  algorithm  in  which  the  next  step  of  propagation  depends  only  on  the  present 
values  of  the  wavefield  in  a  transverse  cross-section  and  the  heterogeneities  between  the  two 
cross-sections.  To  formulate  the  problem,  we  divide  the  crustal  wave  guide  into  a  sequence  of 
vertical  slabs.  The  horizontal  direction  is  chosen  as  the  main  propagation  direction.  The 
geometry  of  the  model  is  shown  in  Figure  3.2a.  Choosing  one  slab  as  the  example.  Figure  3.2b 
shows  the  interaction  between  the  incident  waves  and  the  thin  slab.  By  introducing  the  local 
Bom  approximation,  both  wavefields  and  the  elastic  parameters  can  be  separated  into  two  parts, 
the  background  values  and  the  perturbations.  The  incident  P-  and  S-waves  u£  and  uq  enter  the 
slab  from  the  vertical  plane  atx0 .  After  the  incident  waves  pass  through  the  thin  slab  between  x0 
and  x, ,  interacting  with  the  heterogeneities  within  it,  there  will  be  both  incident  waves  and 
different  types  of  forward  scattered  waves  at  the  exit  plane  at  x, .  The  new  P-wave 
u;  =Uq  +  \JPP  +  Ui7'  is  composed  of  incident  P-wave  and  scattered  P-waves  from  incident  P-  and 
S-waves,  and  the  new  S-wave  us  =  u£  +  U7’5  +  U®  is  composed  of  incident  S-wave  and  scattered 

S-waves  from  incident  P-  and  S-waves.  The  propagation  can  be  solved  using  perturbation  theory. 
The  interaction  between  the  incoming  wave  and  a  heterogeneous  slab  can  be  completed  with  two 
separated  steps.  The  interactions  between  the  incoming  waves  and  the  perturbations  are 
conducted  in  the  spatial  domain  which  give  the  coupling  between  different  wave  types.  Plane 
wave  propagation  through  the  background  medium  is  conducted  in  the  wavenumber  domain.  In 
both  domains,  the  calculations  are  local  and  highly  efficient.  There  is  no  time-consuming  spatial 
or  wavenumber  domain  convolution  involved  The  forward  and  inverse  fast  Fourier  transforms 
alternate  the  wavefield  between  the  two  domains.  By  iteratively  using  this  process  and  making 
the  output  from  one  thin  slab  the  input  of  the  next  thin  slab,  the  wavefield  can  be  propagated 
prograssively  through  the  entire  model. 

The  conventional  wavenumber  integral  method  is  for  horizontally  layered  model  and  the  integral 
is  along  horizontal  wavenumber  kx .  On  the  contrary,  the  elastic  screen  method  for  propagating 
guided  waves  in  crustal  environment  prefers  the  horizontal  as  the  main  propagation  direction. 
Under  this  geometry,  the  postcritical  reflections  result  in  small-angle  dominance  in  wave 
propagation  and  scattering,  which  is  necessary  for  methods  based  on  a  small  angle 
approximation.  Our  velocity  model  is  composed  of  vertical  slabs  and  the  integral  is  along 
vertical  wavenumber  kz .  Wavenumber  integration  has  to  be  modified  for  the  vertical  screen 
geometry. 
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Figure  3.3  Lg  spectral  amplitudes  at  four  different  frequencies  versus  distance  for  crustal  models  with 
topographies  only  (a),  crustal  heterogeneities  only  (b),  and  intrinsic  attenuation  only  (c).  Topography 
with  rms=0.07km  has  a  fluctuation  of  0.87 — 0.75km  in  height.  Qj„  denotes  quality  factor  of  shear 
wave. 


GSP  is  accurate  for  small-angle  propagation  and  scattering  (near  horizontal  for  crustal  wave 
guide  environment).  A  half-space  screen  propagator  has  been  introduced  by  Wu  et  al.  (1996, 
2000)  to  accommodate  the  free-surface  boundary  condition  and  treat  the  SH-wave  propagation  in 
complex  crustal  waveguides.  The  new  one-way  method  for  modeling  regional  SH-waves  has 
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been  calibrated  extensively  with  various  full-wave  methods  for  different  crustal  models,  such  as 
the  wavenumber  integration  method  for  flat  structures  and  finite-difference  method  for 
heterogeneous  crustal  waveguides.  Excellent  agreements  with  these  methods  demonstrated  the 
validity  and  accuracy  of  the  new  one-way  method.  For  a  model  with  propagation  distance  of  250 
km,  dominant  frequency  at  0.5  Hz  and  with  similar  accuracy,  the  GSP  method  is  about  300  times 
faster  than  the  finite-difference  method.  The  method  has  been  applied  to  the  simulation  of  Lg 
propagation  in  random  media  and  the  related  energy  partition  and  attenuation  (Wu  et  al.,  2000). 

It  is  found  that  the  leakage  attenuation  of  Lg  waves  caused  by  forward  large-angle  scattering 
from  random  heterogeneities,  which  scatters  the  guided  waves  out  of  the  trapped  modes  and  into 
the  mantle,  may  contribute  significantly  to  Lg  attenuation  and  blockage  in  some  regions.  The 
apparent  O  for  leakage  attenuation  as  a  function  of  normalized  scale  length  (ka)  of  the  random 
heterogeneities  agrees  well  with  the  scattering  theory.  Later  the  SH  screen  propagator  was 
extended  to  the  case  of  irregular  surface  topography  by  conformal  or  non-conformal  topographic 
transforms  (Wu  et  al.,  1999;  Wu  and  Wu,  2001 ).  In  the  conformal  transform  method,  the 
coordinate  system  is  rotated  according  to  the  local  topographic  slope,  and  the  mirror  image 
method  can  be  applied  to  the  local  plane  surface;  the  non-conformal  method  is  a  surface 
flattening  transform  which  turns  the  free  surface  topography  into  modified  volume  perturbations 
of  elastic  parameters.  The  former  method  is  suitable  to  deal  with  smoothly  varying  topography; 
while  the  non-conformal  method  can  treat  rough  but  moderate  topography. 

In  the  case  of  P-SV  elastic  screen  propagators,  the  derivation  and  application  of  one-way 
screen  propagators  are  much  more  complicated.  Unlike  for  SH-waves,  the  image  method  of 
generating  the  half-space  Green's  function  cannot  be  used  to  account  for  the  effect  of  free 
surface.  Plane  wave  reflection  calculations  are  incorporated  into  the  elastic  screen  method  (Wu  et 
al.,  2000).  Body  waves  including  the  reflected  and  converted  waves  can  be  calculated  by  real 
wavenumber  integration;  while  surface  waves  (Rayleigh  waves)  can  be  done  with  imaginary 
wavenumber  integration. 

3.3  Simulation  of  Lg  Wave  Propagation  and  Blockage  Using  Screen  Method 

Observations  frequently  show  that  regional  Lg  phase  attenuation  is  related  to  many  factors  such 
as  laterally  heterogeneous  velocity  structure,  changes  in  crustal  discontinuities,  crustal  thickness, 
topographic  features,  anelastic  attenuation,  etc.  However,  the  mechanism  of  Lg  wave  attenuation 
has  not  been  fully  understood,  especially  for  tectonically  active  regions  that  have  abnormally  low 
Lg  Q.  For  example  in  the  eastern  Tibetan  Plateau  I  Hz  Lg  Q  is  about  126  (Xie,  2002);  and  it  is 
about  200  -  267  in  various  regions  of  the  western  United  States  (Xie,  1 998),  prior  numerical 
simulations  tend  to  under  estimate  Lg  attenuation  or  blockage.  This  can  be  attributed  to  two 
reasons.  On  one  hand,  realistic  crustal  models  are  quite  complex  so  that  the  parameterized 
models  could  be  poor  approximations.  On  the  other  hand,  the  existing  simulation  approaches 
may  not  have  the  capacity  to  handle  the  effects  of  the  various  factors  mentioned  above. 

In  this  section,  using  the  screen  method  incorporated  with  anelastic  quality  factor  Qin,  Lg  wave 
simulations,  including  the  effects  of  random  topographies,  crustal  heterogeneities  and  anelastic 
attenuation  are  conduced  separately.  We  estimate  apparent  Lg  Q  and  its  power-law  frequency 
dependence  associated  with  different  attenuation  mechanisms.  Then,  we  synthesize  Lg  wave 
responses  to  a  real  crustal  geometry  in  the  Tibetan  Plateau,  assuming  different  heterogeneities 
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and  intrinsic  attenuations.  Preliminary  simulation  results  reveal  high  Lg  attenuation  or  blockage 
after  introducing  low  Qi„  (intrinsic  quality  factor  of  shear  wave)  into  the  model. 


(a)  Frequency  (Hz) 


(b)  Frequency  (Hz) 


Figure  3.4  (a)  Attenuation  coefficient  versus  frequency,  (b)  Apparent  Lg  Q  versus  frequency.  Solid  lines 
correspond  to  the  models  with  random  topographies,  dotted  lines  to  models  with  volume 
heterogeneities,  and  dashed  lines  to  the  anelastic  crustal  waveguides. 

3.3. 1  Surface  and  volume  scattering,  ait  elasticity  and  Lg  O 

The  crustal  model  for  simulations  consists  of  a  crustal  layer  and  a  half-space  mantle.  The 
parameters  of  the  crust  and  mantle  are  =3.5  kmu  ,  pcwi  =  2.*g/cm1,  Vmantlt  =4.5/bw/s,  and 

Pmantu  =3.ig/cm3 .  The  thickness  of  the  crust  is  32  km.  Two  random  models  are  used  for  the  free 
surface  topography.  Both  of  them  have  the  same  correlation  length  of  2.5  km  but  have  different 
RMS  fluctuations.  The  first  has  a  RMS  of  0.07  km  (elevation  between  -0.75  and  0.87  km)  and 
the  second  has  a  RMS  of  0. 14  km  (elevation  between  -1.5  and  1 .75  km).  Random  models  with 
exponential  correlation  function  are  used  for  crust  volume  heterogeneity.  The  correlation  lengths 
are  10  km  in  range  and  5  km  in  depth.  The  RMS  values  used  are  10%  and  20%,  respectively.  A 
source  is  located  at  a  depth  of  8  km  and  has  a  spectrum  function  given  by: 

S(f) - - — T 

1 +(///e)2 

where  /  is  frequency  and  fc  =  2 Hz  is  the  dominant  frequency.  For  the  screen  method,  the  step 
length  for  each  forward  propagation  is  125  m.  The  depth  sampling  interval  is  250  m.  The  number 
of  samples  in  depth  is  5 12.  The  number  of  screens  is  8000  corresponding  to  a  regional  distance 
of  1000km.  We  calculated  6  sets  of  Lg  synthetic  seismograms  for  different  crustal  models. 

Figure  3.3  shows  spectral  amplitudes  versus  distance  for  four  different  frequencies:  0.2,  1.,  2., 
and  4  Hz.  Figure  3.3a  involves  the  scattering  effect  of  irregular  surface  only.  Fluctuation  of  Lg 
spectral  amplitudes  versus  distance  increases  as  frequency  and/or  the  roughness  of  the 
topography  increase,  while  the  average  spectral  amplitudes  show  a  weak  attenuation  for  either 
low  or  high  frequencies.  Figure  3.3b  involves  the  scattering  effect  of  volume  heterogeneities 
only.  Fluctuation  of  Lg  spectral  amplitudes  is  weak  compared  to  the  effect  of  rough  topography. 
Figures  3.3a  and  3.3b  show  that  although  rough  topography  and  volume  heterogeneity  may 
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directly  affect  the  Lg  wave  behavior,  they  cannot  fully  account  for  high  Lg  attenuation  or 
blockage.  Figure  3.3c  involves  only  the  effect  of  anelastic  attenuation.  We  see  that  Lg  spectral 
amplitudes  vary  smoothly  with  distance,  and  high-frequency  Lg  waves  show  strong  attenuation. 


Table  3.1  Estimated  O0  and  r] 


Topography 

Heterogeneity 

Intrinsic  Qjn 

D=0.07km 

D=0. 14km 

B=10% 

B=20% 

500 

300 

Qo 

2029 

945 

1345 

827 

508 

313 

n 

2.320 

1.700 

0.888 

0.73 

0.057 

0.056 

Note:  D  denotes  the  RMS  for  random  topography  and  B  denotes  the  RMS  for 
random  volume  heterogeneity. 


Under  2D  model,  assume  that  Lg  spectral  amplitude  can  be  expressed  as 

Jtf  X 

Al*{f)  =  S(f)e  e’=s(/)C-“ 

where  a  is  the  attenuation  coefficient,  xthe  propagation  distance,  v  the  Lg  group  velocity,  and 
Q  is  the  frequency-dependent  apparent  Lg  Q.  Then  we  can  calculate  attenuation  coefficient  and 
Lg  Q  by  linearly  fitting  Lg  spectral  amplitudes  shown  in  Figure  3.3.  Figure  3.4  shows  the 
estimated  attenuation  coefficients  and  the  frequency  dependence  of  the  apparent  Lg  Q.  Solid 
lines  correspond  to  the  models  with  random  topographies,  dotted  lines  to  the  models  with 
volume  heterogeneities,  and  dashed  lines  to  anelastic  crustal  waveguides.  For  different 
attenuation  mechanisms,  apparent  Lg  Qs  show  different  frequency-dependency. 

The  observed  Lg  Q  usually  fits  a  power-law  frequency  dependence  as: 

0/*(/)=eo/’7 

where  Q0  is  the  Lg  Q  at  1  Hz  reference  frequency  and  n  is  the  power-law  frequency 
dependence.  Using  the  above  relationship  and  Figure  3.4b,  we  can  estimate  Q0  and  n  (listed  in 
Table  1 ).  Table  3. 1  shows  that  the  simulated  Lg  frequency  dependence  tj  for  homogeneous 
anelastic  crusts  is  much  smaller  than  the  observations  from  earthquakes  and  explosions 
( 0.3  <  r]  <  0.6 ),  while  the  j]  associated  with  scattering  by  random  topographies  or  volume 
heterogeneities  varies  strongly.  It  is  possible  to  generate  appropriate  n  by  mixing  the 
contributions  from  both  intrinsic  and  scattering  attenuation  under  the  constraint  of  matching  real 
data. 

3.3.2  Lg  wave  simulation  with  real  crustal  waveguides 

Although  the  propagation  of  Lg  has  been  observed  across  most  of  Asia,  it  has  not  been  observed 
for  paths  crossing  the  Tibetan  Plateau.  It  has  been  demonstrated  that  the  Lg  wave  undergoes 
abnormally  strong  attenuation  in  these  regions.  Figure  3.5a  shows  an  actual  crustal  cross-section 
in  the  Tibetan  Plateau  (Fan  and  Lay,  1998).  This  cross-section  has  severe  changes  in  surface 
topography  and  lateral  velocity  structure.  We  will  introduce  heterogeneity  and  intrinsic 
attenuation  into  the  crustal  model  to  investigate  the  mechanism  of  Lg  attenuation  and  blockage. 
Topography  data  has  a  resolution  of  1  km  in  range  and  Moho  depth  data  has  a  resolution  of  1 0 
km  in  range.  For  purpose  of  numerical  simulation  using  the  screen  method,  both  data  sets  are 
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linearly  interpolated  up  to  a  sampling  interval  of  0. 1 25  km.  The  rest  of  the  parameters  are  the 
same  as  used  in  Figure  3.3.  The  SH  synthetic  seismograms  are  calculated  for  regional  distances 
up  to  1000  km  and  the  Lg  wave  energy  is  obtained  by  summing  up  the  energy  from  all  frequency 
components.  Figure  3.5b  shows  Lg  wave  energy  versus  distance  for  models  with  different 
heterogeneities  and  intrinsic  attenuations.  The  four  lines  from  top  to  bottom  correspond  to:  ( 1 ) 

A  perfect  elastic  homogeneous  waveguide;  (2)  A  perfect  elastic  model  with  10%  RMS  volume 
heterogeneity;  (3)  A  crustal  model  with  10%  RMS  volume  heterogeneity  plus  an  intrinsic  O  of 
500;  and  (4)  A  crustal  model  with  10%  RMS  volume  heterogeneity  plus  an  intrinsic  Q  of  300. 

We  see  that  for  a  perfect  elastic  waveguide,  Lg  energy  fluctuation  closely  follows  the  topography 
and  Lg  energy  decay  is  generally  weak.  Up  to  the  distance  of  1000  km  high-frequency  Lg  phase 
is  still  rather  strong.  After  introducing  intrinsic  attenuation,  Lg  waves  show  higher  attenuation. 


Figure  3.5  (a)  A  cross-section  in  Tibetan  Plateau  (Fan  and  Lay,  1998).  (b)  Lg  nave  energy  versus  distance  for 
crust  models  with  different  small  scale  heterogeneities  and  intrinsic  attenuations,  where  solid  line  is 
for  homogeneous  elastic  waveguide;  dashed  line  is  for  crust  with  10%  RMS  heterogeneity,  dotted  and 
dotted-dashed  lines  are  for  crust  models  with  different  intrinsic  attenuations  (Qjn  =  500  and  300) 
plus  10%  RMS  heterogeneity. 
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Figure  3.6  (a)  Synthetic  seismograms  for  the  crust  with  variable  anelastic  attenuation  along  the  path,  (b) 
synthetic  seismograms  for  the  elastic  crust  with  rms=10%  heterogeneities,  and  (c)  Lg  spectral 
amplitudes  versus  distance  at  4  different  frequencies  calculated  from  (a)  and  (b),  respectively. 

In  the  next  example,  we  use  a  model  with  variable  intrinsic  attenuation  along  the  propagation 
path.  We  assign  the  Qm  to  be  450  for  distances  less  than  300  km;  250  for  distances  between  300 
and  800  km;  and  125  for  distances  greater  than  800  km.  Corresponding  synthetic  seismograms 
are  shown  in  Figure  3.6a.  The  time  axis  is  reduced  by  a  velocity  of  6  km/s.  We  see  that  high- 
frequency  Lg  waves  are  weak  at  800  km  and  absent  at  1000  km.  For  comparison,  synthetic 
seismograms  for  a  perfect  elastic  heterogeneous  waveguide  are  given  in  Figure  3.6b.  High- 
frequency  Lg  waves  are  still  strong  at  distances  up  to  1000  km.  Figure  3.6c  shows  Lg  spectral 
amplitudes  at  four  different  frequencies  versus  distance  calculated  from  synthetic  seismograms 
in  Figures  3.6a  and  b.  The  variable  slope  of  the  high-frequency  Lg  spectral  amplitudes  at 
different  distances  results  from  the  variable  intrinsic  Q  along  the  propagation  path. 
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Figure  3.7  Synthetic  P-SV  Lg  wave  attenuation  curves  for  the  Flora-Asncs  crust  model.  The 
intrinsic  attenuation  in  the  crust  plays  an  important  role  for  Lg  energy  attenuation. 


3.3.3  P-SV  Lg-wave  propagation  in  the  Flora-Asnes  crust  model. 

For  the  P-SV  wave  case,  we  improved  the  stability  of  our  complex  screen  propagator  for  long 
distance  propagation.  At  the  same  time,  we  are  working  on  the  coupling  between  the  body  waves 
and  the  surface  wave.  Figure  3.7  shows  the  P-SV  Lg  wave  attenuation  versus  distance.  The  1-D 
background  model  is  the  Flora-Asnes  crust  model  in  the  NORSAR  region.  A  10%  RMS  random 
velocity  fluctuation  with  correlation  lengths  similar  to  the  SH  wave  case  is  added  to  the 
background  model.  Three  sets  of  P-  and  S-wave  quality  factors  are  used  in  the  calculation  The 
results  show  that  the  intrinsic  attenuation  plays  an  important  role  in  the  Lg  wave  energy  decay. 
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4.  P-SV  BOUNDARY  ELEMENT  METHOD  FOR  REGIONAL  WAVE  PROPAGATION 
4.1  Background 

Regional  phases  have  long  been  recognized  as  important  in  the  study  of  large-scale  crustal 
structures,  small-scale  crustal  heterogeneities,  seismic  sources,  and  determination  of 
underground  explosions  and  earthquakes.  To  improve  path  corrections  as  well  as  verify 
empirical  observations,  numerical  simulations  are  needed  to  estimate  the  path  effects  on  regional 
wave  propagation  A  number  of  regional  numerical  modeling  methods  have  been  developed  to 
model  Lg  propagation  behavior  and  path  effects,  depending  on  the  complexity  of  the  crust 
heterogeneity  in  the  region  considered.  Finite-difference  (FD)  and  finite-element  (FE)  methods 
are  universal  numerical  techniques  for  simulation  in  general  heterogeneous  media.  Because  of 
the  computational  intensity,  these  full-waveform  methods  are  often  limited  to  the  case  of  low 
frequencies  for  Lg  simulation  at  regional  distances.  To  reduce  computation  time,  some 
alternative  and  flexible  approaches,  for  example,  ray  diagrams  (e  g.,  Bostock  and  Kennett,  1990; 
Keers  et  at.,  1996)  and  dynamic-ray-tracing  method  (Gibson  and  Campillo,  1994),  have  been 
used  to  obtain  an  intuitive  understanding  of  path  effects  on  Lg  propagation.  For  large  range  Lg 
wave  propagation,  Wu  et  al.  (1996,  2000)  introduced  a  half-space  GSP  (Generalized  Screen 
Propagator)  for  modeling  the  main  characteristics  of  Lg  (2-D  SH  case)  in  heterogeneous  crustal 
waveguides.  Later  the  method  was  extened  to  include  the  case  of  irregular  surface  topography 
(Wu  and  Wu,  2001) 

Based  on  observations  and  numerical  experiments,  the  large-scale  crustal  structures  with 
variations  in  surface  and  Moho  topography  control  the  principal  characteristics  of  Lg 
propagation.  In  contrast  to  the  FD  and  FE  methods,  the  boundary  integral  methods  are  more 
suitable  for  modeling  complex  reflection  transmission  effects  across  large-scale  crustal  structures 
with  rugged  surface/interface.  For  instance,  the  BE  (boundary  element)  method  provides  a 
geometrically  accurate  description  of  irregular  interfaces.  Because  the  BE  method  is  formulated 
in  terms  of  integrals  along  boundaries,  the  traction-free  condition  for  rugged  free  surface  is  easy 
and  natural  to  treat  in  an  accurate  and  stable  manner.  However,  for  regional  waveguides  up  to  a 
thousand  kilometers,  the  BE  method  in  the  original  form  leads  to  very  large  sized  matrixes  and 
its  inversion  and  implementation  are  prohibitively  expensive.  In  a  previous  paper  (Fu  and  Wu, 
2001),  an  SH  wavefield  connection  technique  is  developed  by  which  the  BE  method  can  be  used 
section  by  section  for  wave  propagation  at  a  far  regional  distance  up  to  several  thousand 
kilometers.  The  method  takes  the  output  of  the  previous  section  as  the  input  of  the  next  section  in 
order  to  complete  the  entire  crustal  waveguide  computation.  The  full-wave  BE  modeling  is 
implemented  within  each  section.  The  wavefield  connection  technique  is  used  to  couple  the 
fields  calculated  in  two  adjoining  sections.  The  division  of  sections  is  based  on  the  complexity  of 
crustal  structures  with  a  criterion  of  minimizing  the  possible  multiple  back  scattering  between 
sections.  This  section-by-section  approach  for  long  regional  waveguides  leads  to  significant 
computational  savings  in  time  and  memory  compared  with  the  whole  waveguide  BE  method. 

The  section-by-section  approach  also  leads  to  a  hybrid  modeling  scheme  of  BE  and  GSP  The 
BE  method  is  implemented  in  the  frequency  domain  and  has  a  kernel  function  compatible  with 
the  GSP  method.  In  the  hybrid  scheme,  the  time-consuming  BE  method  can  be  used  to  handle 
the  sections  with  complicated  boundary  structures  and  severe  surface  topography.  Subsequently 
the  output  will  be  used  as  the  input  to  the  GSP  method  for  modeling  the  sections  with  a  large 
volume  of  moderately  heterogeneous  media  and  mild  topography.  The  hybrid  method  has  been 
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applied  to  two  crustal  waveguide  models  from  the  Tibet  region,  one  with  Lg  blockage  and 
another  without  blockage  (Fu  and  Wu,  2001 ). 

The  object  of  this  section  is  to  develop  the  P-SV  wavefield  connection  technique  for  simulating 
elastic  wave  propagation  in  regional  crustal  waveguides.  We  first  introduce  integral  equation 
representation  for  the  crustal  waveguide  problems.  We  then  briefly  describe  the  principle  of  the 
elastic  BE  method  and  validate  the  computation  programs  using  previously  published  results  We 
develop  P-SV  wavefield  connection  technique  and  test  it  using  numerical  experiments.  Finally, 
we  apply  the  section-by-section  approach  to  Lg  propagation  simulations  in  long  regional 
waveguides. 


Free  Surface 


Moho 


a 


r2  aJ 


Figure  4. 1  Geometry  of  a  simplified  crustal  waveguide  with  irregular  topographic  and  Moho  interface 


4.2  Boundary  Integral  Equation  for  Crustal  Waveguide 

Consider  2-D  steady-state  elastic  wave  propagation  in  a  simplified  crustal  waveguide 
ill  bounded  by  a  rough  free  surface  T i  and  an  irregular  Moho  interface  T2  .  Figure  4. 1  depicts 
the  geometry  of  the  problem.  The  waveguide  medium  is  isotropic  and  homogeneous,  described 
by  the  Lame  constants  and  density  (Xi,  pi,  pi ).  The  displacement  vector  u(r)  at  a  location  r(x,z) 
satisfies  the  following  elastic  wave  equation 

//V‘u(r)  +  (l  +  //)VVu(r)  +  par  u(r)  =  -/?f(r,  to)  (4. 1) 

where  f(r ,(o)  is  the  body  force  occupying  a  region  Qs  u(r)  also  satisfies  the  traction-free 
boundary  condition  on  Ti  and  the  continuity  condition  of  displacement  and  traction  across  the 
Moho.  The  medium  in  the  mantle  Q2  is  described  by  ^ ,  P2 ,  and  P2  .  We  add  two  artificial 
boundaries  at  the  two  truncated  edges  of  the  waveguide  to  form  a  closed  solution  domain 
T  =  T,  +  r2  +  r_ .  Based  on  the  representation  theorem  (e  g.,  Aki  and  Richards,  1980),  equation 
(4.1)  can  be  transformed  into  the  following  integral  formulation: 

CWuC^  +  lHuCrOECr.rO-tCr^GCr.rOk/r^  |npf(r’,<y)G(r,rVr' ,  (4.2) 

where  t(r)  is  the  traction  vector,  the  coefficient  C(r)  generally  depends  on  the  local  geometry  at 
r ,  and  Z(r,r')  and  G(r,r’)  are  the  fundamental  solutions  (Green’s  tensors)  for  displacements  and 
tractions,  respectively.  C(r),£(r,r')and  G(r,  r’)  are  2*2  matrices  for  2-D  problems. 
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where 
¥  = 


Hg\ksr)  -M  H?\ksr)  -fLH\'\k  pr)\ 
ksr\  a  ) 


^H<t\kpr)-HV(ksr) 


ks=a)/p  and  kp  =(o/a,  yp  =(rV-r;.)/r 

In  this  equation,  r  is  the  position  of  the  observation  points  and  r'  is  the  position  of  the  scattering 
points  For  simplicity,  the  source  distribution  consists  of  a  point  source  at  r()  located  inside  fi| 

The  source  integral  over  Qs  in  the  right  side  of  equation  (4.2)  can  be  reduced  to 

\apG{r,r')t{r',a»dr'=  pG(r,r0)f(<y)  (4.3) 

where  f(  ft)  is  the  source  spectrum  vector. 

The  displacement  Green’s  tensor  G(r,  r’)  satisfies: 

[(i  +  //)V  V  •  +//V2  ]G(r,  r' )  +  po)1  G(r,  r’)  =  -8(r  -  r' )  (4.4) 

with  the  solution  given  by  Morse  and  Feshbach  (1953).  The  traction  Green’s  tensor  Z(r,  r')  can 
be  derived  ffomG(r,r’)  using  Hooke’s  law.  The  boundary  integral  representation  for  wave 
propagation  allows  natural  satisfaction  of  Sommerfeld  radiation  boundary  conditions  that  are 
imposed  on  the  far  field  behavior  at  infinity.  No  waves  come  back  to  Q,  through  ,  that  is,  the 
following  integral  on  T,  for  the  interior  problem  vanishes: 

J r_ P-(T> r')u(r') -G(r, r')t(r,r')]dr'=  0,  re  fi,  (4.5) 

For  numerical  calculations,  truncating  the  waveguide  has  to  be  done  and  infinite  boundary 
element  absorbing  boundary  technique  (Fu  and  Wu,  2000)  should  be  applied  to  the  elements  at 
the  truncating  points  on  T,  and  T2 .  Considering  equation  (4.3)  and  (4.5),  and  applying  the 
traction-free  condition  to  the  free  surface  T, ,  equation  (4.2)  for  the  interior  problem  in  the  crust 
(Region  Q, )  is  simplified  to 

C(r)u(r)  +  J  E(r,r')u(r')]dr'+Jr^  L(r,  r')u(r')  -  G(r,  r'  )t(r')]</r’=  ^G(r,r0)f(<w)  (4.6) 

Equation  (4.6)  is  a  starting  point  for  numerical  implementation  for  wave  propagation  simulation. 
In  order  to  solve  u(r)  and  t(r)  on  T2 ,  we  must  build  the  corresponding  boundary  integral 
equation  in  Q2  ■  The  following  integral  formulation  can  be  established  for  reQ2  bounded  by  a 
closed  surface  V  =  T2  +  r_ : 

C(r  )u(r)  +  jr  £(r,  r’  )u(r’ )  -  G  (r,  r’  )t(r'  )]Jr*  =0  (4.7) 

Similarly,  the  integration  over  the  transparent  artificial  boundary  r_  vanishes  for  the  interior 
problem,  simplifying  equation  (4.7)  as 

C(r)u(r)  +  J  E(r, r' )u(r' )  -  G(r, r*  )t(r’  )]cfr*  =  0  (4  8) 
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Equations  (4.7)  and  (4.8)  provide  a  description  of  the  field  through  the  crustal  waveguide, 
making  possible  the  simultaneous  evaluation  of  the  unknowns  (u(r)  on  T i  and  )  u(r)  and  t(r)  on 
T2  )  by  using  the  continuity  of  displacement  and  traction  across  T2  . 

4.3  Boundary  Element  Method  for  Elastic  Wave  Simulation 

The  boundary  integral  representations  described  above  can  be  used  to  calculate  the  field  at  any 
point  inside  £2|  once  u(r)  and  t(r)  are  known  on  the  boundaries.  We  use  the  BE  method  to  solve 
the  values  of  u(r)  and  t(r)  on  the  boundaries.  We  discretize  T i  into  Li  elements  and  N\  nodes, 
and  T2  into  1,2  elements  and  Afe  nodes.  The  total  node  number  is  N.  By  using  linear  interpolation 
shape  functions  <p{%)  in  an  element  between  the  nodes  h  and  /2,  the  variables  u(r)  and  t(r)  are 

approximated  with  the  linear  combination  of  their  node  values  over  the  element,  for  example, 

(see  Fig  4. 1 ) 

u(£)  =  Zu(r'M(£)’  (49> 


where  %  and  /  denote  the  local  coordinate  and  local  node  index  of  an  element.  Using  the 
following  Kronecker  delta  function  notation  relating  the  local  node  code  /  of  an  element  to  the 
global  node  code  j 

[0,  l*j 

/.,  <4I0) 


we  have 


u(#)  =  ZSu(r/M(W 

j- 1  w, 

Letting  C,  =  C(r, )  and  u;  =  u(r, ) ,  equation  (6)  for  /=!  to  N  is  discretized  into 


(4.11) 


“  Z  Z  [f  c<r,  • r'  <#>M  (#)<*•'«)}*,«,  J  -  pC(r„r,  )/<«,) 


which  can  be  further  rewritten  as 

f[h<'>u<»+h>f-gf.<1>]=^G(r„r0)f((B), 

J- 1 

with 


(4.12) 

(4.13) 


u  h 


hT  =££{,.  [2(r„r'(^)M(^Mr-(^)K  +  C«\  (4.14) 

*=1  /-/,  1  * 

h?  =  ££ fr  [E(r„r'(#)M(^r'(#)H  +Cf<5,,  (4.15) 

«=1  (=/,  * 
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(4.16) 


L,  / 


tf '  -  £  £  Jr  [G(r, .  r'  (f  )M  (D*1  <£)]*„ , 


<M  /=/, 


In  this  expression  u1,0  is  the  displacement  vector  on  the  free  surface  F i  and  u*2)  and  t^2)  are  the 

displacement  and  traction  vector  on  the  Moho  r2 .  The  coefficient  matrices  ,  h^2)  and  gj2) , 

obtained  by  numerically  integrating  the  product  of  the  Green’s  tensors  with  interpolation  shape 
functions  over  elements  denotes  a  concentrated  force  generated  at  the  /th  scattering  point  on  T 
and  applied  at  the  /th  observation  point.  For  1  =  /  to  N ,  equation  (4. 13)  can  be  further  compacted 
as  a  matrix  equation: 


H(V‘>  +  H(2)u(2)-G(2)t(2)  =  F,  (4.17) 

where  F  =  />G(r, ,r0)f(6>) .  Similarly,  equation  (4.8)  for  1  =  /  to  N  can  be  discretized  and  compacted 
as 

H(2)u(2)-G(2)t(2)  =0,  (4.18) 

where  the  coefficient  matrices,  H(2)  and  G<2),  are  calculated  using  the  medium  properties  of  the 
mantle  Q2  •  A  simultaneous  system  of  matrix  equations  (4.7)  and  (4.8)  can  be  assembled  using 
the  continuity  of  displacement  and  traction  across  T2  , 

Ju(1)  =u(2) 

I  t(1)  =  t(2) 

by  which  the  unknowns  (u(r)  on  T 1  and  u(r)  and  t(r)  on  T2  )  can  be  solved. 

These  matrices  in  equations  (4  17)  and  (4. 1 8)  are  full  of  complex  coefficients  which  are 
functions  of  frequency,  material  property  and  geometry.  The  BE  method  described  above  can  be 
directly  extended  to  complex  geological  structures  (Fu,  1996)  for  exploration-oriented  seismic 
modeling.  Since  a  large  amount  of  matrix  operations  are  involved  and  the  matrix  for  each 
frequency  component  must  be  inverted,  the  BE  method  is  computationally  intensive  at  high 
frequencies  for  far  regional  waveguides.  To  improve  computation  speed,  a  variable  element 
dimension  technique  can  be  adopted  in  the  program  implementation  (Fu  and  Mu,  1994). 

Bouchon  el  al.  (1995)  suggested  an  approach  of  threshold  criterion  to  make  the  coefficient 
matrices  sparser  by  removing  very  small  entries  in  the  BE  coefficient  matrices.  An  efficient 
modification  to  the  BE  method  can  be  made  using  the  section-by-section  approach  with  the 
wavefield  connection  technique  developed  in  this  study. 

The  computation  program  of  the  elastic  wave  BE  method  described  above  is  tested  by 
dimensionless  frequency  responses  of  a  semi-circular  canyon  with  radius  a.  The  previously 
published  results  (Sanchez-Sesma  et  al.,  1985;  Sanchez-Sesma  and  Campillo,  1991)  for  this 
typical  topographic  structure  are  used  for  comparison.  The  two  sharp  edges  at  x  =  ±a  provide  a 
crucial  target  to  validate  various  numerical  methods.  Figure  4.2  shows  the  comparison  with  good 
agreement  in  both  horizontal  and  vertical  amplitudes  between  our  results  (solid  lines)  and  the 
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Sanchez-Sesma  &  Campillo’s  solutions  for  vertically  incident  plane  P  waves  and  for  various 
normalized  frequencies  T|  defined  as  r\  =  l Ha  ,  where  /  is  the  wavelength.  Poisson’s  ratio  is 
assumed  to  be  1/3.  We  can  see  some  minor  departures  because  of  different  element 
approximations  used  by  these  two  numerical  methods.  Fewer  elements  per  wavelength  will 
reduce  the  size  of  the  resultant  coefficient  matrices.  To  determine  an  applicable  element  number 
per  wavelength,  the  comparisons  for  0°  incidence  and  //  =  1.0  are  given  in  Figure  4.3  for 
different  discretization  rates  of  points  per  wavelength  to  discretize  the  arc  of  the  semicircular 
canyon.  We  see  that  sampling  at  three  points  per  wavelength  could  be  used  for  general 
applications. 


Distance{X/a) 


Figure  4.2:  The  horizontal  and  vertical  amplitude  response  by  our  method  (solid  line)  and  Sanchcz- 
Scsma  &  Campillo  (1991)  (dotted  line)  of  a  semicircular  canyon  topography  to  vertically 
incident  P  waves  for  various  dimensionless  frequency  77 
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Figure  4.3:  The  horizontal  and  vertical  amplitude  responses  of  a  semicircular  canyon  topography  to 
vertically  incident  P  waves  and  TJ  =1.0  with  the  solid  lines  for  21  points  per  wavelength  and 
the  dots  for  various  sampling  rates. 


The  comparisons  above  confirm  the  validity  of  our  formulation  and  computation  code.  We 
now  demonstrate  the  applicability  of  our  program  by  synthesizing  wave  propagation  through  a 
single-layer  crustal  waveguide  with  flat  free  surface.  The  homogeneous  waveguide  shown  in 
Figure  4  4  is  100  km  long  and  32  km  thick,  overlaying  a  flat  mantle  half-space.  The  point 
explosive  source  is  located  on  the  left  boundary  at  2.0-km  depth.  Figure  4.5  shows  the  synthetic 
seismograms  calculated  in  the  frequency  range  of  0-4.5  Hz,  with  receivers  at  1-km  spacing  along 
a  vertical  profile  at  the  distance  of  80  km  from  the  source.  We  see  that  for  elastic  wave 
propagation  this  simple  waveguide  with  both  flat  topography  and  flat  Moho  leads  to  complex 
superposition  of  various  waves,  demonstrating  the  development  of  converted  waves  as  well  as 
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the  formation  of  guided  waves  as  repetitive  reflections  at  both  the  free  surface  and  Moho.  In 
terms  of  propagation  paths  (see  Figure  4.4),  we  can  identify  three  major  systems  of  waves 
generated  in  the  waveguide.  As  shown  in  Fig  4.5,  the  direct  F’-wave  from  the  source  to  receivers 
carries  a  major  part  of  the  source  energy  in  this  flat  boundary  waveguide.  The  P-  wave  incident 
upon  the  free  surface  and  the  Moho  at  oblique  angles  generates  two  systems  of  waves,  one  set  off 
at  the  free  surface  and  another  off  at  the  Moho.  These  two  systems  of  waves  and  their  converted 
waves  bounce  back  and  forth  between  the  free  surface  and  the  Moho.  Since  no  scattering 
mechanism  is  present  in  the  waveguide,  the  constructive  interference  of  repeatedly  reflected 
waves  presents  a  checkerboard-like  pattern  (Jih,  1996)  that  adequately  explains  the  formation  of 
crustal  guided  waves  either  as  multiple  reflections  or  as  higher-modes. 


Figure  4.4:  A  single-layer  flat  crustal  waveguide  with  the  explosive  source  at  2-km  depth.  The  wave 
generated  in  the  waveguide  are  grouped  into  three  systems:  the  first  directly  from  the 
source  to  the  receivers  (solid  line),  the  second  off  at  the  free  surface  (dash  lines)  and  the 
third  off  the  Moho  (dotted  lines). 


Figure  4.5:  Synthetic  seismograms  for  receivers  at  1-km  spacing  along  a  vertical  profde  at  80km  for 
a  point  source.  The  horizontal  (left  panel)  and  vertical  (right  panel)  components  are 
computed  in  the  frequency  range  of  0-3.5hz. 
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Figure  4.6:  The  boundary  connection  technique.  (A)  the  diagrams  show  the  connection  formulation. 
(B)  the  vertical  incident  wavefield  along  the  connecting  boundary  ABB’ ,  (C)  comparison  of 
the  horizontal  component  of  the  wavefield  along  CDD’  calculated  directly  from  the  source 
(thick  lines)  and  the  wavcfield  calculated  using  the  connection  technique  (thin  lines).  (D) 
Comparison  of  the  vertical  component.  The  comparison  show  s  the  validity  of  the  connection 
technique. 


4.4  P-SV  Wavefield  Connection  Technique 

We  aim  to  develop  a  section-by-section  approach  for  simulating  wave  propagation  in  regional 
waveguides  to  reduce  the  computation  cost  of  extremely  large  matrix  operations.  Wavefield 
connection  technique  couples  the  fields  between  two  adjacent  sections.  The  connection 
configuration  is  illustrated  in  Figure  4.6(A).  An  artificial  boundary  is  introduced  as  a  wavefield 
connection  boundary  to  a  crustal  waveguide  consisting  of  an  irregular  free  surface  T i  and  an 
interface  H  .  The  waveguide  is  divided  into  four  subdomains,  Qi,  CI2  in  the  crust,  and  Q, ,  Q;in 
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the  mantle.  The  left  boundaries  of  £2i  and  £2j  ,  and  the  right  boundaries  of  £22  and  £22  are 
assumed  to  extend  to  infinity  at  left  and  right  respectively.  The  fields  in  £2i  and  £2,  ,  can  be 
calculated  from  the  first  section.  Its  output  field  u0(r)  on  the  connection  boundary  r^.  will  be 
used  to  satisfy  the  boundary  condition  when  the  BE  method  is  used  to  calculate  wave 
propagation  in  O2  and  £2: .  The  output  field  is  received  along  the  next  connection  interface  TCD, , 
and  will  be  used  as  the  input  to  the  next  propagation.  With  the  initial  field  u0(r)  known  on  the 
connection  boundary  VAB,  now  we  analyze  wave  scattering  in  £22,  £22  in  order  to  build  boundary 
integral  equations  for  wave  propagation  in  these  sections.  Apparently,  all  the  space  points  of  £2,, 
£22  except  on  r,B.  can  only  receive  the  scattered  field.  Therefore,  the  total  field  u(r)on  the 
boundaries  of  Q2  and  £2,  is  composed  of 


u(r)  = 


K(r)  +  u,(r)  re  T, 


Mr) 


AB' 

re  r^, 


(4.19) 


Applying  the  traction-free  condition  to  the  free  surface  T,  ,  the  boundary  integral  equation  can 


be  built  for  re  £2, and  re  TAB, 


C(r)Uj(r)  +  j  I(r,r>i(r')</rI+f  I(r,r’)[u0(r‘)  +  us(r’)]i/r’+  j  Z(r,r')u5(r’)t/r’ 

Jri+r2  Jr \b  JrCD 

=  L  G(r,r,)tJ(r,)i/r'+f  G(r,r’)[(t0(r')  + tJ(r’)]c/r'+ f  G(r,r')tJ(r')t/r'  (4.20) 

Jr2  as 


The  artificial  boundaries  TAB  and  T^  are  assumed  to  be  transparent,  implying  the  outward 
energy  radiation  across  rAB  and  Tcd  should  always  be  in  the  outward  direction  with  no 
reflection  returning  to  £2,  .  Therefore,  there  is  no  energy  contribution  scattering  from  T,B  and 
TCD  ,  that  is, 

f  [Z(r,r' )u s (r' ) - G(r,r' )t 5 (r' )\ir'  =  0,  re  £22  (4.21) 

AB 

and 

f  [Z(r, r' )u, (r' ) - G(r, r' )t , (r')]c/r*  =  0  re  £2;  (4.22) 

Scattering  from  the  artificial  truncated  points.  A,  B,  C,  and  D,  can  be  handled  using  an  infinite 
element  absorbing  boundary  technique  (Fu  and  Wu,  2000).  Therefore,  equation  (4.20)  is  reduced 
to 


C(r)u,(r)  + j^L(r,r,)uJ(r,)<fr,+j^.jmr,r,)aJ(r,)-G(r,r')ti(r,)]</r' 
=  [a:  l(G(r,  r'  )t0  (!•' )  -  I(r,  r'  )u0  (r'  )]t/r' 


(4.23) 


In  equation  (4 .23),  both  the  initial  displacement  and  traction  on  the  connection  boundary  are 
assumed  to  be  known.  u„(r)can  be  calculated  from  the  previous  section,  incident  traction 

t0(r)can  also  be  calculated  from  the  previous  section.  Alternatively,  we  can  use  the  elastic 

Rayleigh  integrals  (Wu,  1989)  which  contain  either  the  displacement  or  the  traction  field. 
We  can  reduce  the  surface  integral  in  the  right  side  of  equation  (4.23)  to  the  elastic  wave 
Rayleigh  integral  that  only  contains  the  initial  displacement  u0(r) : 
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C(r)us(r)+ f  L(r,r' )us(r')c/r'+f  [I(r,r’ )us(r’ J-Gfrr'jt^r’ )]<*•’  =  -2f  I(r, r')u0(r’)]c/r' 

(4.24) 

where  £(r,r')is  related  to  G(r,r')by  Hooke’s  law: 


Z(r,r')  =  /H(VG(r,r'))  +  //(VG(r,r')  +  G(r,r’)V), 


(4.25) 


where  1  is  the  unit  dyadic  and  G(r,r')V  is  the  transpose  of  VG(r,r')  with  respect  to  the 
corresponding  coordinates.  For  analytic  expressions  of  G  and  £  in  isotropic  homogeneous 
elastic  media,  see  Wu  (1994,  Appendix) 

We  see  that  the  unknowns  in  equation  (4.24)  are  us(r)  on  T|  and  uJ(r)&  ts(r)on  T2  .  In  order 
to  solve  us(r)  &  ts(r) ,  we  must  build  the  corresponding  boundary  integral  equation  in  sub- 
domain  Q2: 

C(r)us(r)+  f  E(r,r,)uJ(r,)</r’+[  Z(r,r')[(u0(r')  +  uJ(r,)]Jr'+  f  Z(r,r’)u0(r')Jr’ 

•M  2 

=  f .G(r,i')t,(i')*'+f  G(r,r’)[(t0(r')  +  ts(r’)]Jr'+  f  G(r,r’)t0(r')Jr'  (4.26) 

W12  DD' 


where  a  sufficient  long  boundary  rBB.is  used  with  its  end  set  as  an  infinite  element.  Similarly, 
since  there  is  no  discontinuity  across  TBB.  and  rDD,  ,  we  can  assume  TBB.  and  rDD,  to  be 
transparent  and  reduce  equation  (4.26)  to 

C(r)Uj(r)+ [  [£(r,r,)uJ(r’)-G(r,r,)tJ(r,)]r/r,=  f  [G(r,r')t0(r')-Z(r,r-)u0(r')]Jr’  (4.27) 

Jr2  JrBB. 

The  integration  term  containing  the  initial  traction  t0(r)  in  the  right  side  of  equation  (4.27)  can 
be  removed  using  the  elastic  wave  Rayleigh  integral  representation,  yielding 

C(r)us  (r)  +  f  [£(r, r’  )u ,  (r’ )  -  G (r, r'  )t s  (r'  )JJr'  =  -2  f  I(r,  r')u0 (r'  )Jr'  (4.28) 

Jr2  jtbb. 

The  continuity  of  displacement  and  traction  across  interface  TBD  is  employed  when  equations 
(4.24)  and  (4.28)  are  combined  to  solve  the  problem.  By  solving  the  joint  boundary  integral 
equations  of  £22and  Q2,  we  can  obtain  the  wavefields  us(r)on  H  ,  Uj(r)  and  ts(r)  on  T2 .  The 
observed  field  along  ^CD’  is  calculated  explicitly  from  the  fields  on  the  boundaries. 


4.5  Validation  for  P-SV  Connection  Technique 

To  test  the  validity  of  the  connection  technique  we  present  a  comparison  between  the  wave 
field  obtained  using  the  BE  method  to  directly  calculate  wave  propagation  from  the  source  to  the 
observation  surface  rcD,  and  the  one  calculated  by  the  connection  scheme  (Fig  4.6).  In  both 
cases  the  source  time  function  are  the  same,  with  a  dominant  frequency  of  1  Hz  First  the 
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intermediate  wave  field  (Fig  4.6  (B))  on  r«,  calculated  from  the  source  is  used  as  the  incident 
field.  Then  the  BE  method  is  used  to  calculate  wave  propagation  from  TAB.  to  rcD, .  The 
dominant  arrivals  for  the  incident  field  at  rAB,  consist  of  direct  P  wave,  pS,  pP,  and  multi¬ 
reflection  between  the  two  layers.  More  multiply  reflected  waves  between  the  free  surface  and 
the  interface  can  be  clearly  seen  in  the  seismograms  at  TCD, .  The  excellent  agreement  between 

the  wavefield  calculated  by  the  two  methods  shown  in  Figure  4.6(C)  and  4.6(D)  confirms  the 
validity  of  the  connection  technique. 
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Figure  4.7  Synthetic  seismogram  for  two-layer  crustal  model  with  random  topography.  The  source 
is  located  at  the  depth  of  2  km.  The  topographic  fluctuation  has  an  exponential  correlation 
function  w  ith  a  RMS  of  0.5  km  and  a  correlation  length  of  5km.  We  can  see  both  for- 
scattering  and  backscattering  in  the  seismogram. 

4.6  Numerical  Examples  and  Applications 

We  use  the  BE  connection  technique  given  above  to  calculate  the  wave  fields  of  several  models 
with  rough  topography.  The  basic  model  is  a  two-layer  crust  over  a  half  space.  The  parameters 
for  this  model  are  listed  in  Table  4. 1 .  The  model  has  two  discontinuities  at  the  depth  of  1 5  km 
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and  37  km.  The  receivers  are  along  the  surface.  An  explosive  source  is  located  at  the  depth  of  2 
km,  and  the  source  time  function  is  a  Ricker  wavelet  with  a  dominant  frequency  of  1  Hz.  Figure 
4.7  shows  the  synthetic  seismogram  for  this  model  with  an  irregular  topography.  The  correlation 
length  of  the  topographic  fluctuation  is  5  km  and  the  RMS  amplitude  is  0.5  km.  Comparing  to 
the  wave  field  (Fig.  4.8)  for  the  model  with  a  flat  surface,  we  can  see  that  the  amplitudes  of 
direct  P  wave  and  Rayleigh  wave  are  relatively  small  due  to  the  scattering  effect  of  the  rough 
topography  Both  fore-scattering  and  back-scattering  waves  can  be  seen  in  Fig  4.7. 


Table  4. 1 :  Crust  velocity  model 


Thickness  (km) 

Vp  (km/s) 

Vs  (km/s) 

p(g/cm3) 

15.00 

6.00 

3.46 

2.80 

22.00 

6.51 

3.76 

3.00 

Infinity 

8.05 

4.65 

3.30 

Figure  4.8  Synthetic  seismogram  for  a  model  with  flat  surface.  Multi-reflected  waves  and  Rayleigh 
waves  can  be  clearly  seen  in  the  figure. 


We  calculated  the  energy  attenuation  curves  for  these  models.  Figure  4.9  shows  the  energy 
attenuation  vs.  distance.  The  one  at  the  bottom  is  the  energy  in  the  whole  seismogram.  The 
energy  is  lost  due  to  the  topography  scattering.  We  calculated  random  topography  with 
exponential  and  Gaussian  distribution  with  0.4  km  RMS  and  1.0  km  RMS.  We  can  see  that 
exponential  topography  is  more  effective  than  the  Gaussian  topography  in  attenuating  the  wave 
energy.  In  addition,  the  larger  the  topographic  fluctuation  is,  the  stronger  the  scattering.  Now  we 
apply  our  connection  technique  to  regional  wave  propagation  simulations  in  a  synthetic  model. 
As  shown  in  Figure  4. 1 0(a),  the  model  is  a  laterally  varying  crustal  model  with  a  600  km  range. 
The  model  is  divided  into  three  200  km  segments.  The  first  segment  is  a  200  km  waveguide  with 
a  Gaussian  hill  which  is  given  by 
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where  *0=62  km,  ho=4  km,  <7=9. 129  km.  The  second  section  has  a  waveguide  necking,  the  height 
of  the  step  is  7  km.  The  third  part  is  a  waveguide  with  random  topography  which  has  a  Gaussian 
correlation  function  with  a  correlation  length  of  5  km  and  the  RMS  height  is  0.6  km.  The  point 
source  is  located  at  8  km  depth.  The  receivers  are  along  both  the  free  surface  and  vertical 
profiles.  The  computation  are  in  the  frequency  range  of  0-2.5  Hz.  We  first  calculate  wave 
propagation  from  source  to  produce  the  incident  wave  field  on  the  first  connection  boundary 
T,  (Fig  4  1 0(b)),  the  multi-reflected  waves,  converted  waves,  head  waves  and  Rayleigh  waves 
can  all  be  clearly  seen  in  the  figure.  Subsequently  the  BE  method  is  used  to  calculate  wave 
propagation  in  the  second  section  and  obtain  the  wave  field  on  the  second  connection  boundary 
T2(Fig  4  10(c)).  Finally  we  use  the  wave  field  on  T2  as  the  incident  field  and  calculate  the 
wavefield  in  the  third  section.  Fig  4  10(d)  shows  the  horizontal  wavefield  on  F, ,  we  can  see 

from  the  scattering  effect  of  the  random  topography,  the  Rayleigh  waves  are  mostly  scattered 
and  some  coda  waves  appear.  We  put  the  wavefield  along  the  surface  obtained  from  the  three 
sections  together  in  Fig  4. 10(e).  Fig  4. 10(f)  is  the  energy  attenuation  along  the  surface.  On  the 
surface  of  the  hill,  the  energy  is  relatively  low  compared  to  the  vicinal  area.  Furthermore,  the 
topography  has  great  effect  on  the  energy  attenuation. 


Figure  4.9:  Energy  attenuation  for  earth  models  with  different  topographies.  On  the  top  is  the 

topography  used  in  our  calculation,  the  thick  solid  line  is  a  Gaussian  random  topographic 
with  a  RMS  of  0.5  km,  thin  solid  line  is  a  Gaussian  random  topography  with  RMS  of  0.25 
km,  the  thick  dash  line  is  exponential  random  topography  with  RMS  of  0.5  km  and  the  thin 
dash  line  is  exponential  random  topography  w  ith  rms  of  0.25  km,  while  the  reference  model 
(flat  surface)  is  show  n  in  dotted  line.  On  the  bottom  are  the  energy  attenuation  curses  for 
the  entire  propagation  distance,  the  lines  arc  the  same  mean  as  those  in  the  top  figure. 
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Figure  4.10  Application  to  a  large  model,  (a)  The  model;  (b)  horizontal  wavefield  on  T, ;  (c) 

horizontal  wavefield  on  T,;  (d)  horizontal  wavefield  on  ;  F,  (c)  horizontal  wavefield  along 
the  surface;  and  (d)  energy  attenuation  along  the  surface. 
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